diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS new file mode 100644 index 0000000..2e093bb --- /dev/null +++ b/.github/CODEOWNERS @@ -0,0 +1,4 @@ +# These owners are the default owners for everything in +# the repo. Unless a later match takes precedence, +* @kousuke-nakano + diff --git a/README.md b/README.md index 9af2653..bcf9ad8 100644 --- a/README.md +++ b/README.md @@ -107,25 +107,16 @@ and type `make html`. The document is generated in `docs/_build/html`. `index.html` is the main page. # Reference -K. Nakano et. al in prepareation (2022). +K. Nakano et. al in prepareation (2023). # How to contribute -Work on the development or on a new branch +Please do not directly push your changes to `main` and `devel` branch. +Please create a pull request on GitHub from your forked repository or a new branch (e.g. `devel-#1`). - git merge development # if you work on a new branch. - git push origin development - -Check the next-version version +# How to check the version info. # Confirm the version number via `setuptools-scm` python -m setuptools_scm e.g., 1.1.4.dev28+gceef293.d20221123 -> = v1.1.4 or v1.1.4-alpha(for pre-release) -Add and push with the new tag - - # Push with tag - git tag # e.g., git tag v1.1.4 # Do not forget "v" before the version number! - git push origin development --tags # or to the new branch - -Send a pull request to the main branch on GitHub. diff --git a/setup.cfg b/setup.cfg index 280af04..cfee5fd 100644 --- a/setup.cfg +++ b/setup.cfg @@ -24,22 +24,22 @@ classifiers = zip_safe = False include_package_data = True packages = find: -python_requires = >="3.7.2" +python_requires = >=3.7.2 install_requires = - matplotlib >= "3.1.1" - numpy >= "1.20.1" - pandas >= "1.2.2" - ase >= "3.21.0" - trexio >= "1.2.0" - pymatgen >= "2021.2.16" - pytest >= "5.2.1" - gitpython >= "3.1.27" - pydriller >= "2.1" - basis-set-exchange >= "0.9" - click >= "8.1.3" - tqdm >= "4.36.1" - setuptools_scm >= "7.0.5" - psutil >= "5.0.0" + matplotlib >= 3.1.1 + numpy >= 1.20.1 + pandas >= 1.2.2 + ase >= 3.21.0 + trexio >= 1.2.0 + pymatgen >= 2021.2.16 + pytest >= 5.2.1 + gitpython >= 3.1.27 + pydriller >= 2.1 + basis-set-exchange >= 0.9 + click >= 8.1.3 + tqdm >= 4.36.1 + setuptools_scm >= 7.0.5 + psutil >= 5.0.0 [options.package_data] * = *.txt, *.rst diff --git a/turbogenius/convertfort10_genius.py b/turbogenius/convertfort10_genius.py index c3bf8dc..13a9bb3 100755 --- a/turbogenius/convertfort10_genius.py +++ b/turbogenius/convertfort10_genius.py @@ -13,6 +13,7 @@ # python modules import os import numpy as np +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -77,13 +78,12 @@ def __init__( self.convertfort10.comment_out(parameter="ny") self.convertfort10.comment_out(parameter="nz") else: - # +- 7.5 bohr from the edges. pos = self.io_fort10.f10structure.positions - Lx = np.max(pos[:, 0]) - np.min(pos[:, 0]) + 3.0 - Ly = np.max(pos[:, 1]) - np.min(pos[:, 1]) + 3.0 - Lz = np.max(pos[:, 2]) - np.min(pos[:, 2]) + 3.0 + Lx = np.max(pos[:, 0]) - np.min(pos[:, 0]) + 6.0 + Ly = np.max(pos[:, 1]) - np.min(pos[:, 1]) + 6.0 + Lz = np.max(pos[:, 2]) - np.min(pos[:, 2]) + 6.0 logger.info( - "Lbox is set to +- 1.5 bohr from the edges of the molecules." + "Lbox is set to +- 3.0 bohr from the edges of the molecules." ) logger.info(f"Lx={Lx}, Ly={Ly}, Lz={Lz}") ax = self.grid_size @@ -150,7 +150,7 @@ def run(self, input_name="convertfort10.input", output_name="out_conv"): flags = self.convertfort10.check_results(output_names=[output_name]) assert all(flags) - def check_results(self, output_names: list = ["out_conv"]) -> bool: + def check_results(self, output_names: Optional[list] = None) -> bool: """ Check the result. @@ -159,6 +159,8 @@ def check_results(self, output_names: list = ["out_conv"]) -> bool: Returns: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_conv"] return self.convertfort10.check_results(output_names=output_names) diff --git a/turbogenius/convertfort10mol_genius.py b/turbogenius/convertfort10mol_genius.py index 75221ae..1d1feff 100755 --- a/turbogenius/convertfort10mol_genius.py +++ b/turbogenius/convertfort10mol_genius.py @@ -13,6 +13,7 @@ # python modules import os import numpy as np +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -40,10 +41,10 @@ class Convertfort10mol_genius(GeniusIO): def __init__( self, - fort10="fort.10_in", - add_random_mo=True, - grid_size=0.10, - additional_mo=0, + fort10: str = "fort.10_in", + add_random_mo: bool = True, + grid_size: float = 0.10, + additional_mo: int = 0, ): self.fort10 = fort10 @@ -178,7 +179,7 @@ def run( flags = self.convertfort10mol.check_results(output_names=[output_name]) assert all(flags) - def check_results(self, output_names: list = ["out_mol"]) -> bool: + def check_results(self, output_names: Optional[list] = None) -> bool: """ Check the result. @@ -187,6 +188,8 @@ def check_results(self, output_names: list = ["out_mol"]) -> bool: Return: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_mol"] return self.convertfort10mol.check_results(output_names=output_names) diff --git a/turbogenius/convertpfaff_genius.py b/turbogenius/convertpfaff_genius.py index 155952d..d842155 100755 --- a/turbogenius/convertpfaff_genius.py +++ b/turbogenius/convertpfaff_genius.py @@ -10,6 +10,7 @@ # python modules import os +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -95,7 +96,7 @@ def run( flags = self.convertpfaff.check_results(output_names=[output_name]) assert all(flags) - def check_results(self, output_names: list = ["out_pfaff"]) -> bool: + def check_results(self, output_names: Optional[list] = None) -> bool: """ Check the result. @@ -104,6 +105,8 @@ def check_results(self, output_names: list = ["out_pfaff"]) -> bool: Returns: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_pfaff"] return self.convertpfaff.check_results(output_names=output_names) diff --git a/turbogenius/correlated_sampling_genius.py b/turbogenius/correlated_sampling_genius.py index 459634e..7a7871c 100755 --- a/turbogenius/correlated_sampling_genius.py +++ b/turbogenius/correlated_sampling_genius.py @@ -12,6 +12,7 @@ # python modules import os +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -52,8 +53,10 @@ def __init__( num_walkers: int = -1, # default -1 -> num of MPI process. maxtime: int = 172800, twist_average: bool = False, - kpoints: list = [1, 1, 1, 0, 0, 0], + kpoints: Optional[list] = None, ): + if kpoints is None: + kpoints = [1, 1, 1, 0, 0, 0] self.in_fort10 = in_fort10 self.corr_fort10 = corr_fort10 @@ -244,8 +247,8 @@ def run( def check_results( self, - vmc_output_names: list = ["out_vmc"], - readforward_output_names: list = ["out_readforward"], + vmc_output_names: Optional[list] = None, + readforward_output_names: Optional[list] = None, ) -> bool: """ Check the result. @@ -256,6 +259,10 @@ def check_results( Returns: bool: True if all the runs were successful, False if an error is detected in the files. """ + if vmc_output_names is None: + vmc_output_names = ["out_vmc"] + if readforward_output_names is None: + readforward_output_names = ["out_readforward"] return self.readforward.check_results( output_names=readforward_output_names ) + self.vmc.check_results(output_names=vmc_output_names) diff --git a/turbogenius/lrdmc_genius.py b/turbogenius/lrdmc_genius.py index b9a9bb4..735e69a 100755 --- a/turbogenius/lrdmc_genius.py +++ b/turbogenius/lrdmc_genius.py @@ -12,6 +12,7 @@ # python modules import os +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -34,6 +35,7 @@ class LRDMC_genius(GeniusIO): fort10 (str): fort.10 WF file lrdmcsteps (int): total number of MCMC steps. alat (float): Lattice space (Bohr) + time_branching: interval between two branching steps. (a.u.) etry (float): Trial Energy (Ha) num_walkers (int): The number of walkers, -1 (default) = the number of MPI processes maxtime (int): Maxtime (sec.) @@ -48,14 +50,17 @@ def __init__( fort10: str = "fort.10", lrdmcsteps: int = 100, alat: float = -0.20, + time_branching: float = 0.10, etry: float = 0.0, num_walkers: int = -1, # default -1 -> num of MPI process. maxtime: int = 172800, twist_average: bool = False, - kpoints: list = [1, 1, 1, 0, 0, 0], + kpoints: Optional[list] = None, force_calc_flag: bool = False, - nonlocalmoves: str = "tmove", # tmove, dla, dlatm + nonlocalmoves: str = "dla", # tmove, dla, dlatm ): + if kpoints is None: + kpoints = [1, 1, 1, 0, 0, 0] self.force_calc_flag = force_calc_flag self.twist_average = twist_average @@ -88,7 +93,9 @@ def __init__( self.lrdmc.set_parameter( parameter="alat", value=alat, namelist="&dmclrdmc" ) - + self.lrdmc.set_parameter( + parameter="tbra", value=time_branching, namelist="&dmclrdmc" + ) typereg, npow = get_nonlocalmoves_setting(nonlocalmoves=nonlocalmoves) self.lrdmc.set_parameter( parameter="typereg", value=typereg, namelist="&dmclrdmc" @@ -257,7 +264,7 @@ def store_result( bin_block: int = 10, warmupblocks: int = 2, correcting_factor: int = 2, - output_names: list = ["out_fn"], + output_names: Optional[list] = None, rerun: bool = False, ) -> None: """ @@ -271,6 +278,8 @@ def store_result( output_names (list): a list of output file names rerun (bool): if true, compute energy and force again even if there are energy and force files. """ + if output_names is None: + output_names = ["out_fn"] self.estimated_time_for_1_generation = ( self.get_estimated_time_for_1_generation(output_names=output_names) ) @@ -305,7 +314,7 @@ def compute_energy_and_forces( ) def get_estimated_time_for_1_generation( - self, output_names: list = ["out_fn"] + self, output_names: Optional[list] = None ) -> float: """ This procedure stores estimated_time_for_1_generation. @@ -316,11 +325,13 @@ def get_estimated_time_for_1_generation( Return: float: estimated_time_for_1_generation. """ + if output_names is None: + output_names = ["out_fn"] return self.lrdmc.get_estimated_time_for_1_generation( output_names=output_names ) - def check_results(self, output_names: list = ["out_fn"]) -> bool: + def check_results(self, output_names: Optional[list] = None) -> bool: """ Check the result. @@ -329,6 +340,8 @@ def check_results(self, output_names: list = ["out_fn"]) -> bool: Return: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_fn"] return self.lrdmc.check_results(output_names=output_names) diff --git a/turbogenius/lrdmc_opt_genius.py b/turbogenius/lrdmc_opt_genius.py index b38a561..01365ea 100755 --- a/turbogenius/lrdmc_opt_genius.py +++ b/turbogenius/lrdmc_opt_genius.py @@ -11,6 +11,7 @@ # python modules import os +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -45,6 +46,7 @@ class LRDMCopt_genius(GeniusIO): learning_rate (float): optimization step size, default values=sr:0.05, lr:0.35 regularization (float): regularization parameter alat (float): Lattice space (Bohr) + time_branching (float): interval between two branching steps (a.u.) etry (float): Trial Energy (Ha) nonlocalmoves (str): Treatment of locality approximation, choose from "tmove", "dla", "dlatm" opt_onebody (bool): flag to optimize onebody Jastrow @@ -61,31 +63,35 @@ class LRDMCopt_genius(GeniusIO): def __init__( self, - fort10="fort.10", - lrdmcoptsteps=100, - steps=10, - bin_block=1, - warmupblocks=1, - num_walkers=-1, # default -1 -> num of MPI process. - maxtime=172800, - optimizer="sr", - learning_rate=0.02, - regularization=0.001, - alat=-0.20, - etry=0.0, - nonlocalmoves="tmove", # tmove, dla, dlatm - opt_onebody=True, - opt_twobody=True, - opt_det_mat=False, - opt_jas_mat=True, - opt_det_basis_exp=False, - opt_jas_basis_exp=False, - opt_det_basis_coeff=False, - opt_jas_basis_coeff=False, - twist_average=False, - kpoints=[1, 1, 1, 0, 0, 0], + fort10: str = "fort.10", + lrdmcoptsteps: int = 100, + steps: int = 10, + bin_block: int = 1, + warmupblocks: int = 1, + num_walkers: int = -1, # default -1 -> num of MPI process. + maxtime: int = 172800, + optimizer: str = "sr", + learning_rate: float = 0.02, + regularization: float = 0.001, + alat: float = -0.20, + time_branching: float = 0.10, + etry: float = 0.0, + nonlocalmoves: str = "dla", # tmove, dla, dlatm + opt_onebody: bool = True, + opt_twobody: bool = True, + opt_det_mat: bool = False, + opt_jas_mat: bool = True, + opt_det_basis_exp: bool = False, + opt_jas_basis_exp: bool = False, + opt_det_basis_coeff: bool = False, + opt_jas_basis_coeff: bool = False, + twist_average: bool = False, + kpoints: Optional[list] = None, ): + if kpoints is None: + kpoints = [1, 1, 1, 0, 0, 0] + self.fort10 = fort10 self.twist_average = twist_average self.kpoints = kpoints @@ -188,6 +194,9 @@ def __init__( self.lrdmcopt.set_parameter( parameter="etry", value=etry, namelist="&dmclrdmc" ) + self.lrdmcopt.set_parameter( + parameter="tbra", value=time_branching, namelist="&dmclrdmc" + ) typereg, npow = get_nonlocalmoves_setting(nonlocalmoves=nonlocalmoves) self.lrdmcopt.set_parameter( @@ -343,13 +352,15 @@ def run( flags = self.lrdmcopt.check_results(output_names=[output_name]) assert all(flags) - def store_result(self, output_names: list = ["out_fn_opt"]) -> None: + def store_result(self, output_names: Optional[list] = None) -> None: """ Store results. energy, energy_error, and estimated_time_for_1_generation are stored in this class. Args: output_names (list): a list of output file names """ + if output_names is None: + output_names = ["out_fn_opt"] self.energy, self.energy_error = self.get_energy( output_names=output_names ) @@ -358,7 +369,9 @@ def store_result(self, output_names: list = ["out_fn_opt"]) -> None: ) def plot_energy_and_devmax( - self, output_names: list = ["out_fn_opt"], interactive: bool = True + self, + output_names: Optional[list] = None, + interactive: bool = True, ): """ plot energy and devmax @@ -367,6 +380,8 @@ def plot_energy_and_devmax( output_names (list): a list of output file names interactive (bool): flag for an interactive plot """ + if output_names is None: + output_names = ["out_fn_opt"] self.lrdmcopt.plot_energy_and_devmax( output_names=output_names, interactive=interactive ) @@ -376,7 +391,7 @@ def average( optwarmupsteps: int = 10, graph_plot: bool = False, input_name: str = "datasfn_opt.input", - output_names: list = ["out_fn_opt"], + output_names: Optional[list] = None, ) -> None: """ Average parameters of fort.10 @@ -387,6 +402,8 @@ def average( output_names (list): a list of output file names graph_plot (bool): Flag for plotting a graph """ + if output_names is None: + output_names = ["out_fn_opt"] flags = self.lrdmcopt.check_results(output_names=output_names) assert all(flags) self.lrdmcopt.average_optimized_parameters( @@ -395,7 +412,7 @@ def average( graph_plot=graph_plot, ) - def get_energy(self, output_names: list = ["out_fn_opt"]) -> list: + def get_energy(self, output_names: Optional[list] = None) -> list: """ return energy list @@ -405,10 +422,12 @@ def get_energy(self, output_names: list = ["out_fn_opt"]) -> list: Return: list: a list of history of energies. """ + if output_names is None: + output_names = ["out_fn_opt"] return self.lrdmcopt.get_energy(output_names=output_names) def get_estimated_time_for_1_generation( - self, output_names: list = ["out_fn_opt"] + self, output_names: Optional[list] = None ) -> float: """ This procedure stores estimated_time_for_1_generation. @@ -419,11 +438,13 @@ def get_estimated_time_for_1_generation( Return: float: estimated_time_for_1_generation. """ + if output_names is None: + output_names = ["out_fn_opt"] return self.lrdmcopt.get_estimated_time_for_1_generation( output_names=output_names ) - def check_results(self, output_names: list = ["out_fn_opt"]) -> bool: + def check_results(self, output_names: Optional[list] = None) -> bool: """ Check the result. @@ -432,6 +453,8 @@ def check_results(self, output_names: list = ["out_fn_opt"]) -> bool: Returns: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_fn_opt"] return self.lrdmcopt.check_results(output_names=output_names) def plot_parameters_history(self, interactive: bool = True) -> None: diff --git a/turbogenius/makefort10_genius.py b/turbogenius/makefort10_genius.py index d3e90bf..93af25a 100755 --- a/turbogenius/makefort10_genius.py +++ b/turbogenius/makefort10_genius.py @@ -14,7 +14,7 @@ import os import numpy as np import glob -from typing import Union +from typing import Union, Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -64,24 +64,31 @@ class Makefort10_genius(GeniusIO): def __init__( self, structure_file: str, - supercell: list = [1, 1, 1], + supercell: Optional[list] = None, det_basis_set: Union[str, list] = "cc-pVQZ", jas_basis_set: Union[str, list] = "cc-pVQZ", det_contracted_flag: bool = True, jas_contracted_flag: bool = True, all_electron_jas_basis_set: bool = True, - pseudo_potential: Union[str, None] = None, + pseudo_potential: Union[str, list, None] = None, det_cut_basis_option: bool = False, jas_cut_basis_option: bool = False, det_exp_to_discard: float = 0.00, jastrow_type: int = -6, complex: bool = False, - phase_up: list = [0.0, 0.0, 0.0], - phase_dn: list = [0.0, 0.0, 0.0], + phase_up: Optional[list] = None, + phase_dn: Optional[list] = None, same_phase_up_dn: bool = False, neldiff: int = 0, ): + if supercell is None: + supercell = [1, 1, 1] + if phase_up is None: + phase_up = [0.0, 0.0, 0.0] + if phase_dn is None: + phase_dn = [0.0, 0.0, 0.0] + self.structure_file = structure_file self.supercell = supercell self.det_basis_set = det_basis_set @@ -546,7 +553,7 @@ def run( """ self.makefort10.run(input_name=input_name, output_name=output_name) - def check_results(self, output_names: list = ["out_make"]) -> bool: + def check_results(self, output_names: list = None) -> bool: """ Check the result. @@ -555,6 +562,8 @@ def check_results(self, output_names: list = ["out_make"]) -> bool: Return: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_make"] return self.makefort10.check_results(output_names=output_names) diff --git a/turbogenius/prep_genius.py b/turbogenius/prep_genius.py index 88fed50..74af09b 100755 --- a/turbogenius/prep_genius.py +++ b/turbogenius/prep_genius.py @@ -13,6 +13,7 @@ # python modules import os import numpy as np +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -51,20 +52,28 @@ class DFT_genius(GeniusIO): def __init__( self, fort10: str = "fort.10", - grid_size: list = [0.1, 0.1, 0.1], - lbox: list = [15.0, 15.0, 15.0], + grid_size: Optional[list] = None, + lbox: Optional[list] = None, smearing: float = 0.0, maxtime: int = 172800, memlarge: bool = False, maxit: int = 50, epsdft: float = 1.0e-5, h_field: float = 0.0, - magnetic_moment_list: list = [], + magnetic_moment_list: Optional[list] = None, xc: str = "lda", # lda or lsda twist_average: bool = False, independent_kpoints: bool = False, - kpoints: list = [1, 1, 1, 0, 0, 0], + kpoints: Optional[list] = None, ): + if grid_size is None: + grid_size = [0.1, 0.1, 0.1] + if lbox is None: + lbox = [15.0, 15.0, 15.0] + if magnetic_moment_list is None: + magnetic_moment_list = [] + if kpoints is None: + kpoints = [1, 1, 1, 0.0, 0] self.fort10 = fort10 self.grid_a, self.grid_b, self.grid_c = grid_size @@ -430,7 +439,7 @@ def run( flags = self.prep.check_results(output_names=[output_name]) assert all(flags) - def check_results(self, output_names: list = ["out_prep"]) -> None: + def check_results(self, output_names: Optional[list] = None) -> None: """ Check the result. @@ -439,6 +448,8 @@ def check_results(self, output_names: list = ["out_prep"]) -> None: Returns: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_prep"] return self.prep.check_results(output_names=output_names) def get_mangetic_moments_3d_array(self): # -> numpy.array(XXX) diff --git a/turbogenius/pyturbo/basis_set.py b/turbogenius/pyturbo/basis_set.py index 7f22004..e3d94fa 100755 --- a/turbogenius/pyturbo/basis_set.py +++ b/turbogenius/pyturbo/basis_set.py @@ -15,7 +15,7 @@ import os import re import numpy as np -from typing import Union +from typing import Union, Optional # logger set from logging import getLogger, StreamHandler, Formatter @@ -52,17 +52,36 @@ class Basis_sets: def __init__( self, - nucleus_index: list = [], - shell_ang_mom: list = [], - shell_ang_mom_turbo_notation: list = [], - shell_factor: list = [], - shell_index: list = [], - exponent: list = [], - coefficient: list = [], - coefficient_imag: list = [], - prim_factor: list = [], + nucleus_index: Optional[list] = None, + shell_ang_mom: Optional[list] = None, + shell_ang_mom_turbo_notation: Optional[list] = None, + shell_factor: Optional[list] = None, + shell_index: Optional[list] = None, + exponent: Optional[list] = None, + coefficient: Optional[list] = None, + coefficient_imag: Optional[list] = None, + prim_factor: Optional[list] = None, ): - + if nucleus_index is None: + nucleus_index = [] + if shell_ang_mom is None: + shell_ang_mom = [] + if shell_ang_mom_turbo_notation is None: + shell_ang_mom_turbo_notation = [] + if shell_factor is None: + shell_factor = [] + if shell_index is None: + shell_index = [] + if exponent is None: + exponent = [] + if coefficient is None: + coefficient = [] + if coefficient_imag is None: + coefficient_imag = [] + if prim_factor is None: + prim_factor = [] + + # variables self.nucleus_index = nucleus_index self.shell_ang_mom = shell_ang_mom self.shell_ang_mom_turbo_notation = shell_ang_mom_turbo_notation @@ -73,15 +92,17 @@ def __init__( self.coefficient_imag = coefficient_imag self.prim_factor = prim_factor - logger.debug(self.nucleus_index) - logger.debug(self.shell_ang_mom) - logger.debug(self.shell_ang_mom_turbo_notation) - logger.debug(self.shell_factor) - logger.debug(self.shell_index) - logger.debug(self.exponent) - logger.debug(self.coefficient) - logger.debug(self.coefficient_imag) - logger.debug(self.prim_factor) + logger.debug(f"nucleus_index={self.nucleus_index}") + logger.debug(f"shell_ang_mom={self.shell_ang_mom}") + logger.debug( + f"shell_ang_mom_turbo_notation={self.shell_ang_mom_turbo_notation}" + ) + logger.debug(f"shell_factor={self.shell_factor}") + logger.debug(f"shell_index={self.shell_index}") + logger.debug(f"exponent={self.exponent}") + logger.debug(f"coefficient={self.coefficient}") + logger.debug(f"coefficient_imag={self.coefficient_imag}") + logger.debug(f"prim_factor={self.prim_factor}") # assertion!! assert len(self.exponent) == len(self.coefficient) @@ -476,7 +497,7 @@ def contracted_to_uncontracted(self): @classmethod def parse_basis_sets_from_gamess_format_files( - cls, files: list = [] + cls, files: Optional[list] = None ): # -> cls """ parse basis sets from GAMESS format files @@ -488,6 +509,8 @@ def parse_basis_sets_from_gamess_format_files( Basis_sets: parsed basis sets """ + if files is None: + files = [] texts = [] for file in files: with open(file, "r") as f: @@ -499,7 +522,7 @@ def parse_basis_sets_from_gamess_format_files( @classmethod def parse_basis_sets_from_eCEPP_format_files( - cls, files: list = [] + cls, files: Optional[list] = None ): # -> cls """ parse basis sets from eCECPP format files @@ -511,6 +534,8 @@ def parse_basis_sets_from_eCEPP_format_files( Basis_sets: parsed basis sets. """ + if files is None: + files = [] texts = [] for file in files: with open(file, "r") as f: @@ -522,7 +547,7 @@ def parse_basis_sets_from_eCEPP_format_files( @classmethod def parse_basis_sets_from_texts( - cls, texts: list = [], format: str = "gamess" + cls, texts: Optional[list] = None, format: str = "gamess" ): # -> cls """ parse basis sets from texts (gamess or eCEPP) @@ -535,6 +560,8 @@ def parse_basis_sets_from_texts( Basis_sets: parsed basis sets. """ + if texts is None: + texts = [] nucleus_index = [] shell_ang_mom = [] shell_ang_mom_turbo_notation = [] @@ -625,13 +652,24 @@ class Basis_set: def __init__( self, - shell_ang_mom=[], - shell_index=[], - exponent_list=[], - coefficient_list=[], - coefficient_imag_list=[], + shell_ang_mom: Optional[list] = None, + shell_index: Optional[list] = None, + exponent_list: Optional[list] = None, + coefficient_list: Optional[list] = None, + coefficient_imag_list: Optional[list] = None, ): - + if shell_ang_mom is None: + shell_ang_mom = [] + if shell_index is None: + shell_index = [] + if exponent_list is None: + exponent_list = [] + if coefficient_list is None: + coefficient_list = [] + if coefficient_imag_list is None: + coefficient_imag_list = [] + + # variables self.shell_ang_mom = shell_ang_mom self.shell_index = shell_index self.exponent_list = exponent_list @@ -969,24 +1007,59 @@ class Det_Basis_sets(Basis_sets): def __init__( self, - nucleus_index: list = [], - shell_ang_mom: list = [], - shell_ang_mom_turbo_notation: list = [], - shell_factor: list = [], - shell_index: list = [], - exponent: list = [], - coefficient: list = [], - coefficient_imag: list = [], - prim_factor: list = [], - hyb_nucleus_index: list = [], - hyb_param_num: list = [], - hyb_shell_ang_mom: list = [], - hyb_shell_ang_mom_turbo_notation: list = [], - hyb_prim_index: list = [], - hyb_coefficient: list = [], - hyb_coefficient_imag: list = [], - number_of_additional_hybrid_orbitals: list = [], + nucleus_index: Optional[list] = None, + shell_ang_mom: Optional[list] = None, + shell_ang_mom_turbo_notation: Optional[list] = None, + shell_factor: Optional[list] = None, + shell_index: Optional[list] = None, + exponent: Optional[list] = None, + coefficient: Optional[list] = None, + coefficient_imag: Optional[list] = None, + prim_factor: Optional[list] = None, + hyb_nucleus_index: Optional[list] = None, + hyb_param_num: Optional[list] = None, + hyb_shell_ang_mom: Optional[list] = None, + hyb_shell_ang_mom_turbo_notation: Optional[list] = None, + hyb_prim_index: Optional[list] = None, + hyb_coefficient: Optional[list] = None, + hyb_coefficient_imag: Optional[list] = None, + number_of_additional_hybrid_orbitals: Optional[list] = None, ): + if nucleus_index is None: + nucleus_index = [] + if shell_ang_mom is None: + shell_ang_mom = [] + if shell_ang_mom_turbo_notation is None: + shell_ang_mom_turbo_notation = [] + if shell_factor is None: + shell_factor = [] + if shell_index is None: + shell_index = [] + if exponent is None: + exponent = [] + if coefficient is None: + coefficient = [] + if coefficient_imag is None: + coefficient_imag = [] + if prim_factor is None: + prim_factor = [] + if hyb_nucleus_index is None: + hyb_nucleus_index = [] + if hyb_nucleus_index is None: + hyb_param_num = [] + if hyb_shell_ang_mom is None: + hyb_shell_ang_mom = [] + if hyb_shell_ang_mom_turbo_notation is None: + hyb_shell_ang_mom_turbo_notation = [] + if hyb_prim_index is None: + hyb_prim_index = [] + if hyb_coefficient is None: + hyb_coefficient = [] + if hyb_coefficient_imag is None: + hyb_coefficient_imag = [] + if number_of_additional_hybrid_orbitals is None: + number_of_additional_hybrid_orbitals = [] + super().__init__( nucleus_index=nucleus_index, shell_ang_mom=shell_ang_mom, @@ -1032,17 +1105,39 @@ class Jas_Basis_sets(Basis_sets): def __init__( self, - nucleus_index=[], - shell_ang_mom=[], - shell_ang_mom_turbo_notation=[], - shell_factor=[], - shell_index=[], - exponent=[], - coefficient=[], - coefficient_imag=[], - prim_factor=[], - number_of_additional_hybrid_orbitals=[], + nucleus_index: Optional[list] = None, + shell_ang_mom: Optional[list] = None, + shell_ang_mom_turbo_notation: Optional[list] = None, + shell_factor: Optional[list] = None, + shell_index: Optional[list] = None, + exponent: Optional[list] = None, + coefficient: Optional[list] = None, + coefficient_imag: Optional[list] = None, + prim_factor: Optional[list] = None, + number_of_additional_hybrid_orbitals: Optional[list] = None, ): + if nucleus_index is None: + nucleus_index = [] + if shell_ang_mom is None: + shell_ang_mom = [] + if shell_ang_mom_turbo_notation is None: + shell_ang_mom_turbo_notation = [] + if shell_factor is None: + shell_factor = [] + if shell_factor is None: + shell_index = [] + if exponent is None: + exponent = [] + if coefficient is None: + coefficient = [] + if coefficient_imag is None: + coefficient_imag = [] + if prim_factor is None: + prim_factor = [] + if number_of_additional_hybrid_orbitals is None: + number_of_additional_hybrid_orbitals = [] + + # variables super().__init__( nucleus_index=nucleus_index, shell_ang_mom=shell_ang_mom, diff --git a/turbogenius/pyturbo/convertfort10.py b/turbogenius/pyturbo/convertfort10.py index 523833a..e7cf56c 100755 --- a/turbogenius/pyturbo/convertfort10.py +++ b/turbogenius/pyturbo/convertfort10.py @@ -15,6 +15,7 @@ # python modules import os import re +from typing import Optional # pyturbo modules from turbogenius.pyturbo.namelist import Namelist @@ -45,10 +46,12 @@ class Convertfort10(FortranIO): def __init__( self, - in_fort10="fort.10_in", - out_fort10="fort.10_out", - namelist=Namelist(), + in_fort10: str = "fort.10_in", + out_fort10: str = "fort.10_out", + namelist: Optional[Namelist] = None, ): + if namelist is None: + namelist = Namelist() """ input values @@ -103,7 +106,7 @@ def run( output_name=output_name, ) - def check_results(self, output_names: list = ["out_conv"]) -> bool: + def check_results(self, output_names: Optional[list] = None) -> bool: """ Check the result. @@ -112,6 +115,8 @@ def check_results(self, output_names: list = ["out_conv"]) -> bool: Returns: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_conv"] flags = [] for output_name in output_names: file_check(output_name) diff --git a/turbogenius/pyturbo/convertfort10mol.py b/turbogenius/pyturbo/convertfort10mol.py index e047085..50f70d0 100755 --- a/turbogenius/pyturbo/convertfort10mol.py +++ b/turbogenius/pyturbo/convertfort10mol.py @@ -15,6 +15,7 @@ # python modules import os import re +from typing import Optional # pyturbo modules from turbogenius.pyturbo.namelist import Namelist @@ -33,10 +34,13 @@ class Convertfort10mol(FortranIO): def __init__( self, - in_fort10="fort.10_in", - namelist=Namelist(), + in_fort10: str = "fort.10_in", + namelist: Optional[Namelist] = None, ): + if namelist is None: + namelist = Namelist() + """ input values """ @@ -57,18 +61,24 @@ def __str__(self): def sanity_check(self): pass - def generate_input(self, input_name="convertfort10mol.input"): + def generate_input(self, input_name: str = "convertfort10mol.input"): self.namelist.write(input_name) logger.info(f"{input_name} has been generated.") - def run(self, input_name="convertfort10mol.input", output_name="out_mol"): + def run( + self, + input_name: str = "convertfort10mol.input", + output_name: str = "out_mol", + ): run( turbo_convertfort10mol_run_command, input_name=input_name, output_name=output_name, ) - def check_results(self, output_names=["out_mol"]): + def check_results(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_mol"] flags = [] for output_name in output_names: file_check(output_name) @@ -86,7 +96,7 @@ def check_results(self, output_names=["out_mol"]): return flags @staticmethod - def read_default_namelist(in_fort10="fort.10_in"): + def read_default_namelist(in_fort10: str = "fort.10_in"): convertfort10mol_default_file = os.path.join( pyturbo_data_dir, "convertfort10mol", "convertfort10mol.input" ) @@ -104,12 +114,12 @@ def read_namelist_from_file(file): return namelist @classmethod - def parse_from_default_namelist(cls, in_fort10="fort.10_in"): + def parse_from_default_namelist(cls, in_fort10: str = "fort.10_in"): namelist = cls.read_default_namelist(in_fort10=in_fort10) return cls(in_fort10=in_fort10, namelist=namelist) @classmethod - def parse_from_file(cls, file, in_fort10="fort.10_in"): + def parse_from_file(cls, file, in_fort10: str = "fort.10_in"): namelist = Namelist.parse_namelist_from_file(file) return cls(in_fort10=in_fort10, namelist=namelist) diff --git a/turbogenius/pyturbo/convertpfaff.py b/turbogenius/pyturbo/convertpfaff.py index 97f3e32..ee41e99 100755 --- a/turbogenius/pyturbo/convertpfaff.py +++ b/turbogenius/pyturbo/convertpfaff.py @@ -12,6 +12,9 @@ """ +# python modules +from typing import Optional + # pyturbo modules from turbogenius.pyturbo.fortranIO import FortranIO from turbogenius.pyturbo.utils.env import turbo_convertfortpfaff_run_command @@ -25,7 +28,9 @@ class Convertpfaff(FortranIO): - def __init__(self, in_fort10="fort.10_in", out_fort10="fort.10_out"): + def __init__( + self, in_fort10: str = "fort.10_in", out_fort10: str = "fort.10_out" + ): """ input values @@ -52,10 +57,10 @@ def generate_input(self): def run( self, - rotate_flag=False, - rotate_angle=0, - scale_mean_field=1000, - output_name="out_pfaff", + rotate_flag: bool = False, + rotate_angle: int = 0, + scale_mean_field: int = 1000, + output_name: str = "out_pfaff", ): if not rotate_flag: @@ -76,7 +81,9 @@ def run( logger.info(f"cmd = {cmd}") run(cmd, input_name=None, output_name=output_name) - def check_results(self, output_names=["out_pfaff"]): + def check_results(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_pfaff"] return [True] @staticmethod @@ -89,12 +96,14 @@ def read_namelist_from_file(): @classmethod def parse_from_default_namelist( - cls, in_fort10="fort.10_in", out_fort10="fort.10_out" + cls, in_fort10: str = "fort.10_in", out_fort10: str = "fort.10_out" ): return cls(in_fort10=in_fort10, out_fort10=out_fort10) @classmethod - def parse_from_file(cls, in_fort10="fort.10_in", out_fort10="fort.10_out"): + def parse_from_file( + cls, in_fort10: str = "fort.10_in", out_fort10: str = "fort.10_out" + ): return cls(in_fort10=in_fort10, out_fort10=out_fort10) diff --git a/turbogenius/pyturbo/io_fort10.py b/turbogenius/pyturbo/io_fort10.py index 4d39957..8758344 100755 --- a/turbogenius/pyturbo/io_fort10.py +++ b/turbogenius/pyturbo/io_fort10.py @@ -155,7 +155,7 @@ class IO_fort10: ) f10jasbasis_sym_end_keyword = "New parameters *$" - def __init__(self, fort10, in_place: bool = True): + def __init__(self, fort10: str = "fort.10", in_place: bool = True): self.fort10 = fort10 self.in_place = in_place @@ -271,7 +271,7 @@ def __init__(self, fort10, in_place: bool = True): # properties!! @property def pp_flag(self): - logger.debug("pseudo_check") + # logger.debug("pseudo_check") atomic_numbers = self.f10structure.atomic_numbers valence_electrons = self.f10structure.valence_electrons diff_list = np.array(atomic_numbers) - np.array(valence_electrons) @@ -500,7 +500,7 @@ def read(self): self.crystal_jastrow_flag = False # logger.debug("check3") - logger.debug("End init of IO_fort10") + # logger.debug("End init of IO_fort10") @property def start_lineno(self): @@ -1009,13 +1009,13 @@ def norm_vec_c(self): @property def atomic_numbers(self): - logger.debug("atomic numbers") + # logger.debug("atomic numbers") self.read() return [i.v for i in self.__atomic_numbers] @property def valence_electrons(self): - logger.debug("valence electrons") + # logger.debug("valence electrons") self.read() return [i.v for i in self.__valence_electrons] @@ -1201,9 +1201,9 @@ def read(self): num_twobody, flag_onebody = return_num_twobody_and_flag_onebody( jastrow_type=self.jastrow_type ) - logger.debug(self.jastrow_type) - logger.debug(num_twobody) - logger.debug(flag_onebody) + # logger.debug(self.jastrow_type) + # logger.debug(num_twobody) + # logger.debug(flag_onebody) if num_twobody == 0: # why here the explicit "if" is needed? # because self.__num_param is not always 0! @@ -1513,7 +1513,7 @@ def det_basis_sets(self): a.v if a.v is not None else 1.0 for a in self.__coefficient ], coefficient_imag=[a.v for a in self.__coefficient_imag], - prim_factor=[1.0] * len(self.__atom_label), + prim_factor=[1.0] * len(self.__exponent), # hybrid hyb_nucleus_index=[a.v - 1 for a in self.__hyb_atom_label], hyb_param_num=[a.v for a in self.__hyb_param_num], @@ -1698,6 +1698,9 @@ def replace_mo_coeff_pure_python( old_mo_coefficient, new_mo_coefficient ): if len(old_mo_coeff) != len(new_mo_coeff): + logger.error( + f"len(old_mo_coeff):{len(old_mo_coeff)} != len(new_mo_coeff):{len(new_mo_coeff)}" + ) raise ValueError line_no_list += [coeff.l for coeff in old_mo_coeff] index_list += [coeff.i for coeff in old_mo_coeff] @@ -1708,7 +1711,7 @@ def replace_mo_coeff_pure_python( logger.debug(f"available memory = {available_memory >> 20} MiB") logger.debug(f"fort10 size = {size_of_fort10 >> 20} MiB") - if available_memory * 0.70 > size_of_fort10: + if available_memory * 0.90 > size_of_fort10: logger.debug("available memory > fort10 size") # version with readlines (more memory, but faster) with open(self.fort10, "r") as f: @@ -1716,92 +1719,64 @@ def replace_mo_coeff_pure_python( # """ straightforward, but established way tqdm_out = TqdmToLogger(logger, level=logging.INFO) + + # new way!! a big loop (line_no_list) is iterated only once!! + line_no_index_list_dict = {line_no: [] for line_no in line_no_list} + for i, line_no in enumerate(line_no_list): + line_no_index_list_dict[line_no].append(i) + for line_no in tqdm( list(set(line_no_list)), file=tqdm_out, miniters=int(len(list(set(line_no_list))) / 10), maxinterval=1.0e5, ): + # old way, very slow!! because + # a big loop (line_no_list) is iterated twice!! + """ line_no_index_list = [ i for i, x in enumerate(line_no_list) if x == line_no ] + """ + # new way + line_no_index_list = line_no_index_list_dict[line_no] + r_index_list = [index_list[i] for i in line_no_index_list] + r_new_mo_coeff_list = [ w_mo_coeff_list[i] for i in line_no_index_list ] + + if len(r_index_list) != len(set(r_index_list)): + logger.error("Duplicated r_index!! It should not happen.") + raise ValueError line = lines[line_no].split() + + # old way (explicit for loop) + """ for r_index, r_new_mo_coeff in zip( r_index_list, r_new_mo_coeff_list ): line[r_index] = r_new_mo_coeff - lines[line_no] = " ".join(list(map(str, line))) + "\n" - # """ - - """ to do: not succesfull yet. it works, but very slow. - lines_new = [] - ind_min = np.min(line_no_list) - ind_max = np.max(line_no_list) - lines_new += lines[:ind_min] - - def replace_line(line_no): - line_no_index_list = [ - i for i, x in enumerate(line_no_list) if x == line_no + """ + # new way (implicit for loop) + line = [ + r_new_mo_coeff_list[r_index_list.index(i)] + if i in r_index_list + else l + for i, l in enumerate(line) ] - r_index_list = [index_list[i] for i in line_no_index_list] - r_new_mo_coeff_list = [ - w_mo_coeff_list[i] for i in line_no_index_list - ] - if len(line_no_index_list) > 0: - line = lines[line_no].split() - for r_index, r_new_mo_coeff in zip( - r_index_list, r_new_mo_coeff_list - ): - line[r_index] = r_new_mo_coeff - return " ".join(list(map(str, line))) + "\n" - else: - return lines[line_no] - - @contextlib.contextmanager - def tqdm_joblib(total: Optional[int] = None, **kwargs): - - pbar = tqdm(total=total, miniters=1, smoothing=0, **kwargs) - - class TqdmBatchCompletionCallback( - joblib.parallel.BatchCompletionCallBack - ): - def __call__(self, *args, **kwargs): - pbar.update(n=self.batch_size) - return super().__call__(*args, **kwargs) - - old_batch_callback = joblib.parallel.BatchCompletionCallBack - joblib.parallel.BatchCompletionCallBack = ( - TqdmBatchCompletionCallback - ) - - try: - yield pbar - finally: - joblib.parallel.BatchCompletionCallBack = ( - old_batch_callback - ) - pbar.close() + lines[line_no] = " ".join(list(map(str, line))) + "\n" - logger.warning("mutliprocessing using joblib!!") - with tqdm_joblib(len(range(ind_min, ind_max))): - lines_new += joblib.Parallel(n_jobs=-1)( - joblib.delayed(replace_line)(i) - for i in range(ind_min, ind_max) - ) - lines_new += lines[ind_max:] - lines = lines_new - """ + # """ with open(self.fort10, "w") as f: f.writelines(lines) else: logger.debug("available memory < fort10 size") - # version with readline (less memory, but slower) + + # version with readline (less memory, but super slow!!) sed_fort10 = self.fort10 + "_sed_tmp" with open(self.fort10, "r") as f: with open(sed_fort10, "w") as fw: diff --git a/turbogenius/pyturbo/lrdmc.py b/turbogenius/pyturbo/lrdmc.py index 98d0bfd..5c45258 100755 --- a/turbogenius/pyturbo/lrdmc.py +++ b/turbogenius/pyturbo/lrdmc.py @@ -18,6 +18,7 @@ import os import re import numpy as np +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -46,8 +47,13 @@ class LRDMC(FortranIO): def __init__( - self, in_fort10="fort.10", namelist=Namelist(), twist_average=False + self, + in_fort10: str = "fort.10", + namelist: Optional[Namelist] = None, + twist_average: bool = False, ): + if namelist is None: + namelist = Namelist() """ input values @@ -61,6 +67,36 @@ def __init__( self.namelist = namelist self.twist_average = twist_average + # manual k-grid! [[[kx, ky, kz, wkp for up], ....], [[# kx, ky, kz, wkp for dn], ...]] + self.__manual_kpoints = [] + + @property + def manual_kpoints(self): + return self.__manual_kpoints + + @manual_kpoints.setter + def manual_kpoints(self, kpoints): + assert len(kpoints) == 2 + kpoints_up, kpoints_dn = kpoints + assert len(kpoints_up) == len(kpoints_dn) + for kup, kdn in zip(kpoints_up, kpoints_dn): + assert len(kup) == 4 # kx, ky, kz, wkp for up + assert len(kdn) == 4 # kx, ky, kz, wkp for dn + self.__manual_kpoints = kpoints + self.namelist.set_parameter( + parameter="yes_kpoints", value=".true.", namelist="¶meters" + ) + # self.namelist.set_parameter(parameter="yeswrite10", value=".true.", namelist="&optimization") + self.namelist.set_parameter( + parameter="kp_type", value=2, namelist="&kpoints" + ) + self.namelist.set_parameter( + parameter="nk1", value=len(kpoints_up), namelist="&kpoints" + ) + self.namelist.set_parameter( + parameter="double_kpgrid", value=".true.", namelist="&kpoints" + ) + def __str__(self): output = [ @@ -71,8 +107,42 @@ def __str__(self): def sanity_check(self): pass - def generate_input(self, input_name="datasfn.input"): + def generate_input(self, input_name: str = "datasfn.input"): self.namelist.write(input_name) + # check if twist_average is manual + if self.twist_average == 2: # k-points are set manually + kpoints_up, kpoints_dn = self.manual_kpoints + # read + with open(input_name, "r") as f: + lines = f.readlines() + kpoint_index = [ + True if re.match(r".*&kpoints.*", line) else False + for line in lines + ].index(True) + insert_lineno = -1 + for i, line in enumerate(lines[kpoint_index + 1 :]): + if re.match(r".*/.*", line): + insert_lineno = kpoint_index + 1 + i + 1 + break + if re.match(r".*&.*", line): + insert_lineno = kpoint_index + 1 + i + break + if insert_lineno == -1: + logger.error("No suitable position for KPOINT is found") + raise ValueError + # write dn + lines.insert(insert_lineno, "\n") + for kx, ky, kz, wk in reversed(kpoints_dn): + lines.insert(insert_lineno, f"{kx} {ky} {kz} {wk}\n") + lines.insert(insert_lineno, "\n") + # write up + for kx, ky, kz, wk in reversed(kpoints_up): + lines.insert(insert_lineno, f"{kx} {ky} {kz} {wk}\n") + lines.insert(insert_lineno, "KPOINTS\n") + lines.insert(insert_lineno, "\n") + # saved + with open(input_name, "w") as f: + f.writelines(lines) logger.info(f"{input_name} has been generated.") def run(self, input_name="datasfn.input", output_name="out_fn"): @@ -84,7 +154,9 @@ def run(self, input_name="datasfn.input", output_name="out_fn"): output_name=output_name, ) - def check_results(self, output_names=["out_fn"]): + def check_results(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_fn"] flags = [] for output_name in output_names: file_check(output_name) @@ -98,7 +170,11 @@ def check_results(self, output_names=["out_fn"]): flags.append(False) return flags - def get_estimated_time_for_1_generation(self, output_names=["out_fn"]): + def get_estimated_time_for_1_generation( + self, output_names: Optional[list] = None + ): + if output_names is None: + output_names = ["out_fn"] out_min = [] for output_name in output_names: @@ -124,7 +200,12 @@ def get_estimated_time_for_1_generation(self, output_names=["out_fn"]): return ave_time_1_generation # sec. def get_energy( - self, init=10, correct=10, bin=10, num_proc=-1, rerun=False + self, + init: int = 10, + correct: int = 10, + bin: int = 10, + num_proc: int = -1, + rerun: bool = False, ): force_compute_flag = False if rerun: @@ -152,7 +233,7 @@ def get_energy( return energy, error @staticmethod - def read_energy(twist_average=False): + def read_energy(twist_average: bool = False): if twist_average: line = get_line_from_file(file="pip0_fn.d", line_no=0).split() energy = float(line[3]) @@ -164,7 +245,12 @@ def read_energy(twist_average=False): return energy, error def get_forces( - self, init=10, correct=10, bin=10, num_proc=-1, rerun=False + self, + init: int = 10, + correct: int = 10, + bin: int = 10, + num_proc: int = -1, + rerun: bool = False, ): force_compute_flag = False @@ -275,7 +361,12 @@ def get_forces( return force_matrix, force_matrix_error_bar def compute_energy_and_forces( - self, init=10, correct=10, bin=10, pulay=1, num_proc=-1 + self, + init: int = 10, + correct: int = 10, + bin: int = 10, + pulay: int = 1, + num_proc: int = -1, ): if self.twist_average: if num_proc == -1: @@ -288,7 +379,7 @@ def compute_energy_and_forces( num_proc = os.cpu_count() if num_proc > 1: command = turbo_forcefn_kpoints_para_run_command - #command = turbo_forcefn_kpoints_run_command # for the time being!!! because paperoga does not have sufficient memory. + # command = turbo_forcefn_kpoints_run_command # for the time being!!! because paperoga does not have sufficient memory. else: command = turbo_forcefn_kpoints_run_command cmd = "{:s} {:d} {:d} {:d} {:d} {:d}".format( @@ -310,13 +401,13 @@ def read_default_namelist(): return namelist @staticmethod - def read_namelist_from_file(file): + def read_namelist_from_file(file: str): namelist = Namelist.parse_namelist_from_file(file) return namelist @classmethod def parse_from_default_namelist( - cls, in_fort10="fort.10", twist_average=False + cls, in_fort10: str = "fort.10", twist_average: bool = False ): namelist = cls.read_default_namelist() return cls( @@ -324,7 +415,9 @@ def parse_from_default_namelist( ) @classmethod - def parse_from_file(cls, file, in_fort10="fort.10", twist_average=False): + def parse_from_file( + cls, file: str, in_fort10: str = "fort.10", twist_average: bool = False + ): namelist = Namelist.parse_namelist_from_file(file) logger.debug("parse from file, lrdmc.py") return cls( diff --git a/turbogenius/pyturbo/lrdmcopt.py b/turbogenius/pyturbo/lrdmcopt.py index 46745bb..fccaf2f 100755 --- a/turbogenius/pyturbo/lrdmcopt.py +++ b/turbogenius/pyturbo/lrdmcopt.py @@ -21,6 +21,7 @@ import pandas as pd import matplotlib.pyplot as plt import matplotlib.ticker as ticker +from typing import Optional # turbo-genius modules from turbogenius.pyturbo.namelist import Namelist @@ -47,11 +48,13 @@ class LRDMCopt(FortranIO): def __init__( self, - in_fort10="fort.10", - namelist=Namelist(), - pp_flag=False, - twist_average=False, + in_fort10: str = "fort.10", + namelist: Optional[Namelist] = None, + pp_flag: bool = False, + twist_average: bool = False, ): + if namelist is None: + namelist = Namelist() """ input values @@ -75,11 +78,15 @@ def __str__(self): def sanity_check(self): pass - def generate_input(self, input_name): + def generate_input(self, input_name: str): self.namelist.write(input_name) logger.info(f"{input_name} has been generated. \n") - def run(self, input_name="datasfn_opt.input", output_name="out_fn_opt"): + def run( + self, + input_name: str = "datasfn_opt.input", + output_name: str = "out_fn_opt", + ): run( turbo_qmc_run_command, input_name=input_name, @@ -87,7 +94,9 @@ def run(self, input_name="datasfn_opt.input", output_name="out_fn_opt"): ) remove_file(file="stop.dat") - def check_results(self, output_names=["out_fn_opt"]): + def check_results(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_fn_opt"] flags = [] for output_name in output_names: file_check(output_name) @@ -103,8 +112,10 @@ def check_results(self, output_names=["out_fn_opt"]): return flags def plot_energy_and_devmax( - self, output_names=["out_fn_opt"], interactive=False + self, output_names: Optional[list] = None, interactive: bool = False ): + if output_names is None: + output_names = ["out_fn_opt"] # plot the energies and devmax logger.info("Plotting the energies and devmax") @@ -193,7 +204,9 @@ def plot_energy_and_devmax( ) plt.close() - def get_energy(self, output_names=["out_fn_opt"]): + def get_energy(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_fn_opt"] out_min = [] for output_name in output_names: @@ -223,8 +236,9 @@ def get_energy(self, output_names=["out_fn_opt"]): return energy_list, error_list - def get_devmax(self, output_names=["out_fn_opt"]): - + def get_devmax(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_fn_opt"] out_min = [] for output_name in output_names: with open(output_name, "r") as f: @@ -243,7 +257,11 @@ def get_devmax(self, output_names=["out_fn_opt"]): return devmax_list - def get_estimated_time_for_1_generation(self, output_names=["out_fn_opt"]): + def get_estimated_time_for_1_generation( + self, output_names: Optional[list] = None + ): + if output_names is None: + output_names = ["out_fn_opt"] out_min = [] for output_name in output_names: @@ -270,9 +288,9 @@ def get_estimated_time_for_1_generation(self, output_names=["out_fn_opt"]): def average_optimized_parameters( self, - equil_steps=10, - input_file_used="datasfnopt.input", - graph_plot=False, + equil_steps: int = 10, + input_file_used: str = "datasfnopt.input", + graph_plot: bool = False, ): if self.twist_average: raise NotImplementedError @@ -421,13 +439,13 @@ def read_default_namelist(): return namelist @staticmethod - def read_namelist_from_file(file): + def read_namelist_from_file(file: str): namelist = Namelist.parse_namelist_from_file(file) return namelist @classmethod def parse_from_default_namelist( - cls, in_fort10="fort.10", twist_average=False + cls, in_fort10: str = "fort.10", twist_average: bool = False ): namelist = cls.read_default_namelist() return cls( @@ -435,7 +453,9 @@ def parse_from_default_namelist( ) @classmethod - def parse_from_file(cls, file, in_fort10="fort.10", twist_average=False): + def parse_from_file( + cls, file, in_fort10: str = "fort.10", twist_average: bool = False + ): namelist = Namelist.parse_namelist_from_file(file) return cls( in_fort10=in_fort10, namelist=namelist, twist_average=twist_average diff --git a/turbogenius/pyturbo/makefort10.py b/turbogenius/pyturbo/makefort10.py index 0644619..3c80e1d 100755 --- a/turbogenius/pyturbo/makefort10.py +++ b/turbogenius/pyturbo/makefort10.py @@ -18,6 +18,7 @@ import re import numpy as np from decimal import Decimal +from typing import Optional # set logger from logging import getLogger, StreamHandler, Formatter @@ -42,12 +43,22 @@ class Makefort10(FortranIO): def __init__( self, - structure=Structure(), - det_basis_sets=Det_Basis_sets(), - jas_basis_sets=Jas_Basis_sets(), - pseudo_potentials=Pseudopotentials(), - namelist=Namelist(), + structure: Optional[Structure] = None, + det_basis_sets: Optional[Det_Basis_sets] = None, + jas_basis_sets: Optional[Jas_Basis_sets] = None, + pseudo_potentials: Optional[Pseudopotentials] = None, + namelist: Optional[Namelist] = None, ): + if structure is None: + structure = Structure() + if det_basis_sets is None: + det_basis_sets = Det_Basis_sets() + if jas_basis_sets is None: + jas_basis_sets = Jas_Basis_sets() + if pseudo_potentials is None: + pseudo_potentials = Pseudopotentials() + if namelist is None: + namelist = Namelist() """ input values @@ -68,7 +79,9 @@ def __str__(self): def sanity_check(self): assert self.structure.natom == self.det_basis_sets.nuclei_num - def generate_input(self, input_name, basis_sets_unique_element=True): + def generate_input( + self, input_name: str, basis_sets_unique_element: bool = True + ): # pseudo potential generation (pseudo.dat) self.pseudo_potentials.write_pseudopotential_turborvb_file() @@ -181,7 +194,7 @@ def generate_input(self, input_name, basis_sets_unique_element=True): else: # label = "ATOM_{:f}".format(fake_atomic_number) # :f format does not work. label = "ATOM_{:f}".format( - Decimal(str(fake_atomic_number)).normalize() + Decimal(str(round(fake_atomic_number, 8))).normalize() ) nshelldet = len( [ @@ -390,14 +403,20 @@ def generate_input(self, input_name, basis_sets_unique_element=True): with open(input_name, "a") as f: f.writelines(output) - def run(self, input_name="makefort10.input", output_name="out_make"): + def run( + self, + input_name: str = "makefort10.input", + output_name: str = "out_make", + ): run( turbo_makefort10_run_command, input_name=input_name, output_name=output_name, ) - def check_results(self, output_names=["out_make"]): + def check_results(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_make"] flags = [] for output_name in output_names: file_check(output_name) @@ -412,12 +431,14 @@ def check_results(self, output_names=["out_make"]): @staticmethod def read_default_namelist( - structure=Structure(), - jastrow_type=-6, - neldiff=0, - nel=-1, - complex_flag=False, + structure: Optional[Structure] = None, + jastrow_type: int = -6, + neldiff: int = 0, + nel: int = -1, + complex_flag: bool = False, ): + if structure is None: + structure = Structure() makefort10_default_file = os.path.join( pyturbo_data_dir, "makefort10", "makefort10.input" ) diff --git a/turbogenius/pyturbo/namelist.py b/turbogenius/pyturbo/namelist.py index bf91d71..7c0d7d6 100755 --- a/turbogenius/pyturbo/namelist.py +++ b/turbogenius/pyturbo/namelist.py @@ -15,6 +15,7 @@ # python modules import re +from typing import Optional, Union from collections import OrderedDict # turbo-genius modules @@ -27,9 +28,11 @@ class Namelist: - def __init__(self, namelist=OrderedDict()): - assert type(namelist) == type( - OrderedDict() + def __init__(self, namelist: Optional[OrderedDict] = None): + if namelist is None: + namelist = OrderedDict() + assert isinstance( + namelist, OrderedDict ), "Please use OrderedDict() for the namelist!" self.__namelist = namelist @@ -37,7 +40,12 @@ def __init__(self, namelist=OrderedDict()): def parameters(self): return self.__namelist - def set_parameter(self, parameter, value, namelist=None): + def set_parameter( + self, + parameter: str, + value: Union[int, float, str], + namelist: Optional[dict] = None, + ): if namelist is None: for key, parameters in self.__namelist.items(): if parameter in parameters.keys(): diff --git a/turbogenius/pyturbo/prep.py b/turbogenius/pyturbo/prep.py index 7f37fa5..26d57e9 100755 --- a/turbogenius/pyturbo/prep.py +++ b/turbogenius/pyturbo/prep.py @@ -17,6 +17,7 @@ import os import re import numpy as np +from typing import Optional # pyturbo modules from turbogenius.pyturbo.namelist import Namelist @@ -48,13 +49,23 @@ class Prep(FortranIO): def __init__( self, - in_fort10="fort.10", - namelist=Namelist(), - nelocc_list=[], - neloccdn_list=[], - magnetic_moments_3d_array=[], # dim = 3, shape = (nzs,nys,nxs) - twist_average=False, # False or 0: single-k, True or 1: Monkhorst-Pack, 2: manual k-grid + in_fort10: str = "fort.10", + namelist: Optional[Namelist] = None, + nelocc_list: Optional[list] = None, + neloccdn_list: Optional[list] = None, + magnetic_moments_3d_array: Optional[ + list + ] = None, # dim = 3, shape = (nzs,nys,nxs) + twist_average: bool = False, # False or 0: single-k, True or 1: Monkhorst-Pack, 2: manual k-grid ): + if namelist is None: + namelist = Namelist() + if nelocc_list is None: + nelocc_list = [] + if neloccdn_list is None: + neloccdn_list = [] + if magnetic_moments_3d_array is None: + magnetic_moments_3d_array = [] """ input values @@ -129,7 +140,7 @@ def __str__(self): def sanity_check(self): pass - def generate_input(self, input_name): + def generate_input(self, input_name: str = "prep.input"): """ Args: input_name (str): Name of input file @@ -211,7 +222,9 @@ def run(self, input_name="prep.input", output_name="out_prep"): output_name=output_name, ) - def check_results(self, output_names=["out_prep"]): + def check_results(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_prep"] flags = [] for output_name in output_names: file_check(output_name) @@ -232,7 +245,7 @@ def check_results(self, output_names=["out_prep"]): return flags @staticmethod - def read_default_namelist(in_fort10="fort.10"): + def read_default_namelist(in_fort10: str = "fort.10"): prep_default_file = os.path.join( pyturbo_data_dir, "prep", "prep.input" ) @@ -283,10 +296,12 @@ def read_namelist_from_file(file): @classmethod def parse_from_default_namelist( cls, - in_fort10="fort.10", - magnetic_moments_3d_array=[], - twist_average=False, + in_fort10: str = "fort.10", + magnetic_moments_3d_array: Optional[list] = None, + twist_average: bool = False, ): + if magnetic_moments_3d_array is None: + magnetic_moments_3d_array = [] namelist = cls.read_default_namelist(in_fort10=in_fort10) # fort10 io_fort10 = IO_fort10(fort10=in_fort10) @@ -320,7 +335,7 @@ def parse_from_default_namelist( ) @classmethod - def parse_from_file(cls, file, in_fort10="fort.10"): + def parse_from_file(cls, file: str, in_fort10: str = "fort.10"): namelist = Namelist.parse_namelist_from_file(file) logger.warning( f"nelocc_list and neloccdn_list are not read from {file}" diff --git a/turbogenius/pyturbo/pseudopotentials.py b/turbogenius/pyturbo/pseudopotentials.py index bef55c4..76c3ae1 100755 --- a/turbogenius/pyturbo/pseudopotentials.py +++ b/turbogenius/pyturbo/pseudopotentials.py @@ -19,6 +19,7 @@ import shutil import random import string +from typing import Optional # turbo-genius modules from turbogenius.pyturbo.utils.env import pyturbo_tmp_dir, turborvb_bin_root @@ -39,17 +40,36 @@ class Pseudopotentials: def __init__( self, - max_ang_mom_plus_1=[], - z_core=[], - cutoff=[], - nucleus_index=[], - element_list=[], - ang_mom=[], - exponent=[], - coefficient=[], - power=[], + max_ang_mom_plus_1: Optional[list] = None, + z_core: Optional[list] = None, + cutoff: Optional[list] = None, + nucleus_index: Optional[list] = None, + element_list: Optional[list] = None, + ang_mom: Optional[list] = None, + exponent: Optional[list] = None, + coefficient: Optional[list] = None, + power: Optional[list] = None, ): - + if max_ang_mom_plus_1 is None: + max_ang_mom_plus_1 = [] + if z_core is None: + z_core = [] + if cutoff is None: + cutoff = [] + if nucleus_index is None: + nucleus_index = [] + if element_list is None: + element_list = [] + if ang_mom is None: + ang_mom = [] + if exponent is None: + exponent = [] + if coefficient is None: + coefficient = [] + if power is None: + power = [] + + # variables self.element_list = element_list self.max_ang_mom_plus_1 = max_ang_mom_plus_1 self.z_core = z_core @@ -60,15 +80,15 @@ def __init__( self.coefficient = coefficient self.power = power - logger.debug(self.element_list) - logger.debug(self.max_ang_mom_plus_1) - logger.debug(self.z_core) - logger.debug(self.cutoff) - logger.debug(self.nucleus_index) - logger.debug(self.ang_mom) - logger.debug(self.exponent) - logger.debug(self.coefficient) - logger.debug(self.power) + logger.debug(f"element_list={self.element_list}") + logger.debug(f"max_ang_mom_plus_1={self.max_ang_mom_plus_1}") + logger.debug(f"z_core={self.z_core}") + logger.debug(f"cutoff={self.cutoff}") + logger.debug(f"nucleus_index={self.nucleus_index}") + logger.debug(f"ang_mom={self.ang_mom}") + logger.debug(f"exponent={self.exponent}") + logger.debug(f"coefficient={self.coefficient}") + logger.debug(f"power={self.power}") # assertion!! assert len(set(self.nucleus_index)) == len(self.max_ang_mom_plus_1) @@ -86,7 +106,7 @@ def nuclei_num(self): def ecp_num(self): return len(self.ang_mom) - def write_pseudopotential_turborvb_file(self, file="pseudo.dat"): + def write_pseudopotential_turborvb_file(self, file: str = "pseudo.dat"): with open(file, "w") as f: output = [] output.append("ECP\n") @@ -149,7 +169,7 @@ def write_pseudopotential_turborvb_file(self, file="pseudo.dat"): f.writelines(output) - def set_cutoffs(self, tollerance=0.00001): + def set_cutoffs(self, tollerance: float = 0.00001): logger.debug(self.cutoff) new_cutoff = [] current_dir = os.getcwd() @@ -244,7 +264,9 @@ def set_cutoffs(self, tollerance=0.00001): os.chdir(current_dir) @classmethod - def parse_pseudopotential_from_turborvb_pseudo_dat(cls, file="pseudo.dat"): + def parse_pseudopotential_from_turborvb_pseudo_dat( + cls, file: str = "pseudo.dat" + ): max_ang_mom_plus_1 = [] z_core = [] @@ -408,15 +430,23 @@ def parse_pseudopotential_from_turborvb_format_files(cls, files): class Pseudopotential: def __init__( self, - max_ang_mom_plus_1=0, - z_core=0, - element=None, - cutoff=0.0, - ang_mom=[], - exponent=[], - coefficient=[], - power=[], + max_ang_mom_plus_1: int = 0, + z_core: int = 0, + element: Optional[str] = None, + cutoff: float = 0.0, + ang_mom: Optional[list] = None, + exponent: Optional[list] = None, + coefficient: Optional[list] = None, + power: Optional[list] = None, ): + if ang_mom is None: + ang_mom = [] + if exponent is None: + exponent = [] + if coefficient is None: + coefficient = [] + if power is None: + power = [] self.element = element self.max_ang_mom_plus_1 = max_ang_mom_plus_1 @@ -445,14 +475,14 @@ def ecp_num(self): return len(self.ang_mom) @classmethod - def parse_pseudopotential_from_gamess_format_file(cls, file): + def parse_pseudopotential_from_gamess_format_file(cls, file: str): with open(file, "r") as f: text = f.readlines() text = "".join(text) return cls.parse_pseudopotential_from_gamess_format_text(text=text) @classmethod - def parse_pseudopotential_from_gamess_format_text(cls, text): + def parse_pseudopotential_from_gamess_format_text(cls, text: str): """ " http://myweb.liu.edu/~nmatsuna/gamess/input/ECP.html -card 1- PNAME, PTYPE, IZCORE, LMAX+1 @@ -537,7 +567,7 @@ def parse_pseudopotential_from_gamess_format_text(cls, text): ) @classmethod - def parse_pseudopotential_from_turborvb_format_file(cls, file): + def parse_pseudopotential_from_turborvb_format_file(cls, file: str): with open(file, "r") as f: text = f.readlines() text = "".join(text) @@ -551,7 +581,7 @@ def parse_pseudopotential_from_turborvb_format_file(cls, file): @classmethod def parse_pseudopotential_from_turborvb_format_test( - cls, text, z_core, element + cls, text: str, z_core: int, element: str ): ang_mom = [] exponent_list = [] @@ -598,14 +628,14 @@ def parse_pseudopotential_from_turborvb_format_test( ) @classmethod - def parse_pseudopotential_from_eCEPP_format_file(cls, file): + def parse_pseudopotential_from_eCEPP_format_file(cls, file: str): with open(file, "r") as f: text = f.readlines() text = "".join(text) return cls.parse_pseudopotential_from_eCEPP_format_text(text=text) @classmethod - def parse_pseudopotential_from_eCEPP_format_text(cls, text): + def parse_pseudopotential_from_eCEPP_format_text(cls, text: str): ang_mom = [] exponent_list = [] coefficient_list = [] diff --git a/turbogenius/pyturbo/readforward.py b/turbogenius/pyturbo/readforward.py index a688187..55e814b 100755 --- a/turbogenius/pyturbo/readforward.py +++ b/turbogenius/pyturbo/readforward.py @@ -16,6 +16,7 @@ # python modules import os import re +from typing import Optional # pyturbo modules from turbogenius.pyturbo.namelist import Namelist @@ -35,9 +36,11 @@ class Readforward(FortranIO): def __init__( self, - in_fort10="fort.10", - namelist=Namelist(), + in_fort10: str = "fort.10", + namelist: Optional[Namelist] = None, ): + if namelist is None: + namelist = Namelist() """ input values @@ -58,18 +61,24 @@ def __str__(self): def sanity_check(self): pass - def generate_input(self, input_name="readforward.input"): + def generate_input(self, input_name: str = "readforward.input"): self.namelist.write(input_name) logger.info(f"{input_name} has been generated.") - def run(self, input_name="datasvmc.input", output_name="out_readforward"): + def run( + self, + input_name: str = "datasvmc.input", + output_name: str = "out_readforward", + ): run( turbo_readforward_run_command, input_name=input_name, output_name=output_name, ) - def check_results(self, output_names=["out_readforward"]): + def check_results(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_readforward"] flags = [] for output_name in output_names: file_check(output_name) @@ -87,7 +96,7 @@ def check_results(self, output_names=["out_readforward"]): return flags @staticmethod - def read_default_namelist(in_fort10="fort.10"): + def read_default_namelist(in_fort10: str = "fort.10"): readforward_default_file = os.path.join( pyturbo_data_dir, "readforward", "readforward.input" ) @@ -95,17 +104,17 @@ def read_default_namelist(in_fort10="fort.10"): return namelist @staticmethod - def read_namelist_from_file(file): + def read_namelist_from_file(file: str): namelist = Namelist.parse_namelist_from_file(file) return namelist @classmethod - def parse_from_default_namelist(cls, in_fort10="fort.10"): + def parse_from_default_namelist(cls, in_fort10: str = "fort.10"): namelist = cls.read_default_namelist(in_fort10=in_fort10) return cls(in_fort10=in_fort10, namelist=namelist) @classmethod - def parse_from_file(cls, file, in_fort10="fort.10"): + def parse_from_file(cls, file, in_fort10: str = "fort.10"): namelist = Namelist.parse_namelist_from_file(file) return cls(in_fort10=in_fort10, namelist=namelist) diff --git a/turbogenius/pyturbo/structure.py b/turbogenius/pyturbo/structure.py index 24d634b..a805bd3 100755 --- a/turbogenius/pyturbo/structure.py +++ b/turbogenius/pyturbo/structure.py @@ -17,6 +17,7 @@ import math import numpy as np from numpy import linalg as LA +from typing import Optional # python material modules from ase import Atoms @@ -36,15 +37,21 @@ class Cell: def __init__( self, - vec_a=np.array([0.0, 0.0, 0.0], dtype=float), # bohr - vec_b=np.array([0.0, 0.0, 0.0], dtype=float), # bohr - vec_c=np.array([0.0, 0.0, 0.0], dtype=float), # bohr + vec_a: Optional[list] = None, + vec_b: Optional[list] = None, + vec_c: Optional[list] = None, ): + if vec_a is None: + vec_a = [0.0, 0.0, 0.0] # bohr + if vec_b is None: + vec_b = [0.0, 0.0, 0.0] # bohr + if vec_c is None: + vec_c = [0.0, 0.0, 0.0] # bohr # cell vectors (the units are bohr) - self.__vec_a = vec_a - self.__vec_b = vec_b - self.__vec_c = vec_c + self.__vec_a = np.array(vec_a, dtype=float) + self.__vec_b = np.array(vec_b, dtype=float) + self.__vec_c = np.array(vec_c, dtype=float) # calc norm self.__norm_vec_a = LA.norm(self.__vec_a) @@ -232,11 +239,19 @@ def parse_cell_from_celldm( class Structure: def __init__( self, - cell=Cell(), - atomic_numbers=[], - element_symbols=[], - positions=np.array([[]]), # 3 * N matrix, the unit is bohr!! + cell: Optional[Cell] = None, + atomic_numbers: Optional[list] = None, + element_symbols: Optional[list] = None, + positions: Optional[list] = None, # 3 * N matrix, the unit is bohr!! ): + if cell is None: + cell = Cell() + if atomic_numbers is None: + atomic_numbers = [] + if element_symbols is None: + element_symbols = [] + if positions is None: + positions = np.array([[]]) self.__cell = cell self.__atomic_numbers = atomic_numbers diff --git a/turbogenius/pyturbo/utils/downloader.py b/turbogenius/pyturbo/utils/downloader.py index 9faa12e..8a36007 100755 --- a/turbogenius/pyturbo/utils/downloader.py +++ b/turbogenius/pyturbo/utils/downloader.py @@ -10,108 +10,125 @@ """ -#python modules +# python modules import os -import sys import re -import numpy as np -import urllib from urllib.request import urlopen -import xml.etree.ElementTree as ET -import re -from pathlib import Path -from ase.data import chemical_symbols -from time import sleep import git import tempfile import pathlib -import pydriller import itertools import time -import random from tqdm import tqdm from tqdm.contrib import itertools as titertools -from collections import OrderedDict +from typing import Optional -#pymatgen +# pymatgen from pymatgen.core.periodic_table import Element -#basissetexchange +# ASE +from ase.data import chemical_symbols + +# basissetexchange import basis_set_exchange as bse -#turbo-genius modules -sys.path.append(os.path.dirname(os.path.abspath(__file__))) -from env import pyturbo_data_dir, pyturbo_root +# turbo-genius modules +from turbogenius.pyturbo.utils.env import pyturbo_root + +# Logger +from logging import getLogger, StreamHandler, Formatter + +logger = getLogger("pyturbo").getChild(__name__) -#Logger -from logging import config, getLogger, StreamHandler, Formatter -logger = getLogger('pyturbo').getChild(__name__) class ccECP: C_URL = "https://github.com/QMCPACK/pseudopotentiallibrary.git" list_of_basis_all = [f"cc-pV{s}Z" for s in "DTQ56"] list_of_basis_all += [f"ang-{x}" for x in list_of_basis_all] - def __init__(self, basis_sets_output_dir=None, pseudo_potential_output_dir=None): + def __init__( + self, + basis_sets_output_dir: Optional[str] = None, + pseudo_potential_output_dir: Optional[str] = None, + ): if basis_sets_output_dir is None: - self.basis_sets_output_dir=None + self.basis_sets_output_dir = None else: - self.basis_sets_output_dir=basis_sets_output_dir + self.basis_sets_output_dir = basis_sets_output_dir if pseudo_potential_output_dir is None: self.ecp_output_dir = None else: self.ecp_output_dir = pseudo_potential_output_dir - def to_file(self, element_list=None, basis_list=None): + def to_file( + self, + element_list: Optional[list] = None, + basis_list: Optional[list] = None, + ): + if element_list is None: + element_list = [] + if basis_list is None: + basis_list = [] if self.basis_sets_output_dir is not None: os.makedirs(self.basis_sets_output_dir, exist_ok=True) if self.ecp_output_dir is not None: os.makedirs(self.ecp_output_dir, exist_ok=True) elements = {} - def add_data(p, tempdir, elements = elements): - if self.basis_sets_output_dir is None: return + + def add_data(p, tempdir, elements=elements): + if self.basis_sets_output_dir is None: + return element = str(p.parent.parent.name) typ = str(p.parent.name) - typ_core = typ.replace('ccECP', '') + typ_core = typ.replace("ccECP", "") logger.info(typ) logger.info(typ_core) - element_path = p.parent.parent - + # element_path = p.parent.parent + typ = str(p.parent.name) - typ_path = str(p.parent.name) + # typ_path = str(p.parent.name) if element_list is not None: if element not in element_list: return - #Load Basis sets - basis = [] + # Load Basis sets + # basis = [] for r in (p.parent).glob("**/*"): if re.match("[A-z]{1,2}\.[A-z\-]*cc-.*\.gamess", r.name): if re.match("[A-z]{1,2}\.[A-z\-]*cc-.*\.gamess", r.name): - basis_name = re.match("[A-z]{1,2}\.([A-z\-]*cc-.*)\.gamess", r.name).group(1) + basis_name = re.match( + "[A-z]{1,2}\.([A-z\-]*cc-.*)\.gamess", r.name + ).group(1) if basis_list is not None: if basis_name not in basis_list: continue logger.info(element) name = f"{basis_name}{typ_core}.basis" with open(r, "r") as fhandle: - with open(os.path.join(self.basis_sets_output_dir, f"{element}_{name}"), "w") as fhandle_out: + with open( + os.path.join( + self.basis_sets_output_dir, + f"{element}_{name}", + ), + "w", + ) as fhandle_out: fhandle_out.write(fhandle.read()) - + tempdir = pathlib.Path(tempfile.mkdtemp()) git.Repo.clone_from(self.C_URL, tempdir) - #logger.debug(list((tempdir/"recipes").glob("**/*"))) - for p in (tempdir/"recipes").glob("**/*"): + # logger.debug(list((tempdir/"recipes").glob("**/*"))) + for p in (tempdir / "recipes").glob("**/*"): if re.match("[A-z]{1,2}\.ccECP\.gamess", p.name): logger.debug(p) add_data(p, tempdir) - if self.ecp_output_dir is None: return + if self.ecp_output_dir is None: + return element = str(p.parent.parent.name) typ = str(p.parent.name) - typ_core=typ.replace('ccECP','') + typ_core = typ.replace("ccECP", "") logger.info(typ) logger.info(typ_core) @@ -119,71 +136,122 @@ def add_data(p, tempdir, elements = elements): if element not in element_list: continue logger.info(p.name) - logger.info(p.name.replace('.ccECP.gamess', "")) + logger.info(p.name.replace(".ccECP.gamess", "")) with open(p, "r") as fhandle: - with open(os.path.join(self.ecp_output_dir, f"{p.name.replace('.ccECP.gamess', '_ccECP'+typ_core+'.pseudo')}"), "w") as fhandle_out: + with open( + os.path.join( + self.ecp_output_dir, + f"{p.name.replace('.ccECP.gamess', '_ccECP'+typ_core+'.pseudo')}", + ), + "w", + ) as fhandle_out: fhandle_out.write(fhandle.read()) - def all_to_file(self, sleep_time=1): - self.to_file(element_list=chemical_symbols, basis_list=self.list_of_basis_all) + def all_to_file(self, sleep_time: float = 1): + self.to_file( + element_list=chemical_symbols, basis_list=self.list_of_basis_all + ) + class BSE: - #URL: https://www.basissetexchange.org + # URL: https://www.basissetexchange.org list_of_basis_all = [f"cc-pV{s}Z" for s in "DTQ56"] list_of_basis_all += [f"ang-{x}" for x in list_of_basis_all] - def __init__(self, basis_sets_output_dir=None, pseudo_potential_output_dir=None): + def __init__( + self, + basis_sets_output_dir: Optional[str] = None, + pseudo_potential_output_dir: Optional[str] = None, + ): if basis_sets_output_dir is None: - self.basis_sets_output_dir=None + self.basis_sets_output_dir = None else: - self.basis_sets_output_dir=basis_sets_output_dir - - def to_file(self, element_list, basis_list, sleep_time=1.0): - if self.basis_sets_output_dir is None: return + self.basis_sets_output_dir = basis_sets_output_dir + + def to_file( + self, + element_list: Optional[list] = None, + basis_list: Optional[list] = None, + sleep_time: int = 1.0, + ): + if element_list is None: + element_list = [] + if basis_list is None: + basis_list = [] + if self.basis_sets_output_dir is None: + return os.makedirs(self.basis_sets_output_dir, exist_ok=True) for e, b in itertools.product(element_list, basis_list): - #time.sleep(sleep_time + random.randint(sleep_time)) + # time.sleep(sleep_time + random.randint(sleep_time)) try: - basis_text=bse.get_basis(b, elements=e, fmt='gamess_us') - E=Element(e) - pat = re.compile("^.*?DATA.*" + E.long_name + "\n(.+)\$END", re.DOTALL | re.IGNORECASE) + basis_text = bse.get_basis(b, elements=e, fmt="gamess_us") + E = Element(e) + pat = re.compile( + "^.*?DATA.*" + E.long_name + "\n(.+)\$END", + re.DOTALL | re.IGNORECASE, + ) m = pat.match(basis_text) - if (m): - bas=m.group(1) - with open(os.path.join(self.basis_sets_output_dir, f"{e}_{b}.basis"), "w") as fhandle: + if m: + bas = m.group(1) + with open( + os.path.join( + self.basis_sets_output_dir, f"{e}_{b}.basis" + ), + "w", + ) as fhandle: fhandle.write(bas) except KeyError: - logger.debug(f"element={e}, basis={b} do not exist in the database.") + logger.debug( + f"element={e}, basis={b} do not exist in the database." + ) - def all_to_file(self, sleep_time=1): - self.to_file(element_list=chemical_symbols, basis_list=self.list_of_basis_all, sleep_time=sleep_time) + def all_to_file(self, sleep_time: float = 1): + self.to_file( + element_list=chemical_symbols, + basis_list=self.list_of_basis_all, + sleep_time=sleep_time, + ) class BFD: U = "http://burkatzki.com/pseudos/step4.2.php?format=gamess&element={e}&basis={b}" list_of_basis_all = [f"v{s}z" for s in "dtq56"] - #list_of_basis_all += [f"{x}_ano" for x in list_of_basis_all] + # list_of_basis_all += [f"{x}_ano" for x in list_of_basis_all] - def __init__(self, basis_sets_output_dir=None, pseudo_potential_output_dir=None): + def __init__( + self, + basis_sets_output_dir: Optional[str] = None, + pseudo_potential_output_dir: Optional[str] = None, + ): if basis_sets_output_dir is None: - self.basis_sets_output_dir=None + self.basis_sets_output_dir = None else: - self.basis_sets_output_dir=basis_sets_output_dir + self.basis_sets_output_dir = basis_sets_output_dir if pseudo_potential_output_dir is None: self.ecp_output_dir = None else: self.ecp_output_dir = pseudo_potential_output_dir - def to_file(self, element_list, basis_list, sleep_time=1.5): + def to_file( + self, + element_list: Optional[list] = None, + basis_list: Optional[list] = None, + sleep_time: float = 1.5, + ): + if element_list is None: + element_list = [] + if basis_list is None: + basis_list = [] if self.basis_sets_output_dir is None and self.ecp_output_dir is None: return if self.basis_sets_output_dir is not None: os.makedirs(self.basis_sets_output_dir, exist_ok=True) if self.ecp_output_dir is not None: os.makedirs(self.ecp_output_dir, exist_ok=True) - if "X" in element_list: element_list.remove("X") + if "X" in element_list: + element_list.remove("X") with tqdm(titertools.product(element_list, basis_list)) as pbar: - for i, (e,b) in enumerate(pbar): + for i, (e, b) in enumerate(pbar): logger.debug(f"b={b:s}") logger.debug(f"[BFD] e={e:s} b={b:s}") @@ -191,21 +259,47 @@ def to_file(self, element_list, basis_list, sleep_time=1.5): if self.basis_sets_output_dir is not None: if self.ecp_output_dir is not None: flag_downloaded = ( - (os.path.isfile(os.path.join(self.ecp_output_dir, f"{e}_BFD.pseudo")) or - os.path.isfile(os.path.join(self.ecp_output_dir, f"{e}_BFD.NaN"))) and - (os.path.isfile(os.path.join(self.basis_sets_output_dir, f"{e}_{b}.basis")) or - os.path.isfile(os.path.join(self.basis_sets_output_dir, f"{e}_{b}.NaN"))) + os.path.isfile( + os.path.join( + self.ecp_output_dir, f"{e}_BFD.pseudo" + ) + ) + or os.path.isfile( + os.path.join( + self.ecp_output_dir, f"{e}_BFD.NaN" + ) + ) + ) and ( + os.path.isfile( + os.path.join( + self.basis_sets_output_dir, + f"{e}_{b}.basis", + ) + ) + or os.path.isfile( + os.path.join( + self.basis_sets_output_dir, f"{e}_{b}.NaN" + ) + ) ) else: - flag_downloaded = ( - (os.path.isfile(os.path.join(self.basis_sets_output_dir, f"{e}_{b}.basis")) or - os.path.isfile(os.path.join(self.basis_sets_output_dir, f"{e}_{b}.NaN"))) + flag_downloaded = os.path.isfile( + os.path.join( + self.basis_sets_output_dir, f"{e}_{b}.basis" + ) + ) or os.path.isfile( + os.path.join( + self.basis_sets_output_dir, f"{e}_{b}.NaN" + ) ) else: if self.ecp_output_dir is not None: - flag_downloaded = ( - (os.path.isfile(os.path.join(self.ecp_output_dir, f"{e}_BFD.pseudo")) or - os.path.isfile(os.path.join(self.ecp_output_dir, f"{e}_BFD.NaN"))) + flag_downloaded = os.path.isfile( + os.path.join( + self.ecp_output_dir, f"{e}_BFD.pseudo" + ) + ) or os.path.isfile( + os.path.join(self.ecp_output_dir, f"{e}_BFD.NaN") ) else: @@ -215,56 +309,97 @@ def to_file(self, element_list, basis_list, sleep_time=1.5): logger.debug(f"b={b:s}") pbar.set_description(f"[BFD] e={e:s} b={b:s}") time.sleep(sleep_time) - logger.debug(self.U.format(b = b, e = e)) - source=urlopen(self.U.format(b = b, e = e)).read() + logger.debug(self.U.format(b=b, e=e)) + source = urlopen(self.U.format(b=b, e=e)).read() source = str(source) logger.debug(source) - pat=re.compile("^.*?(" + e + "\s0.*$)",re.M|re.DOTALL) + pat = re.compile("^.*?(" + e + "\s0.*$)", re.M | re.DOTALL) x = pat.sub("\g<1>", source) x = re.sub("\", "\n", x) x = re.sub("\ ", "", x) x = re.sub("\ ", "", x) x = re.sub(".*html.*$", "", x) - pat=re.compile("^.*?(" + e + "\s\d.*)("+e+".*QMC.*)$",re.M|re.DOTALL) + pat = re.compile( + "^.*?(" + e + "\s\d.*)(" + e + ".*QMC.*)$", + re.M | re.DOTALL, + ) m = pat.match(x) - if(m): + if m: bas = m.group(1) - bas="\n".join([i for i in bas.split("\n")[1:]]).rstrip("\n") # remove the first line. - ecp = m.group(2) # ECP is OK. + bas = "\n".join( + [i for i in bas.split("\n")[1:]] + ).rstrip( + "\n" + ) # remove the first line. + ecp = m.group(2) # ECP is OK. if self.ecp_output_dir is not None: - with open(os.path.join(self.ecp_output_dir, f"{e}_BFD.pseudo"), "w") as fo: - #logger.info(f"ecp -> {os.path.join(self.ecp_output_dir, f'{e}_BFD.pseudo')}") + with open( + os.path.join( + self.ecp_output_dir, f"{e}_BFD.pseudo" + ), + "w", + ) as fo: + # logger.info(f"ecp -> {os.path.join(self.ecp_output_dir, f'{e}_BFD.pseudo')}") fo.write(ecp) if self.basis_sets_output_dir is not None: if len(bas) > 15: - with open(os.path.join(self.basis_sets_output_dir, f"{e}_{b}.basis"), "w") as fo: - #logger.info(f"basis -> {os.path.join(self.basis_sets_output_dir, f'{e}_{b}.basis')}") + with open( + os.path.join( + self.basis_sets_output_dir, + f"{e}_{b}.basis", + ), + "w", + ) as fo: + # logger.info(f"basis -> {os.path.join(self.basis_sets_output_dir, f'{e}_{b}.basis')}") fo.write(bas) else: - with open(os.path.join(self.basis_sets_output_dir, f"{e}_{b}.NaN"), "w") as fo: + with open( + os.path.join( + self.basis_sets_output_dir, + f"{e}_{b}.NaN", + ), + "w", + ) as fo: fo.write("NA in database") else: - with open(os.path.join(self.ecp_output_dir, f"{e}_BFD.NaN"), "w") as fo: + with open( + os.path.join(self.ecp_output_dir, f"{e}_BFD.NaN"), + "w", + ) as fo: fo.write("NA in database") - with open(os.path.join(self.basis_sets_output_dir, f"{e}_{b}.NaN"), "w") as fo: + with open( + os.path.join( + self.basis_sets_output_dir, f"{e}_{b}.NaN" + ), + "w", + ) as fo: fo.write("NA in database") else: - logger.info(f"e={e:s} b={b:s} already downloaded or NaN in the database") + logger.info( + f"e={e:s} b={b:s} already downloaded or NaN in the database" + ) - def all_to_file(self, sleep_time=1.5): + def all_to_file(self, sleep_time: float = 1.5): logger.debug(self.list_of_basis_all) - self. to_file(basis_list=self.list_of_basis_all, element_list=chemical_symbols, sleep_time=sleep_time) + self.to_file( + basis_list=self.list_of_basis_all, + element_list=chemical_symbols, + sleep_time=sleep_time, + ) + if __name__ == "__main__": logger = getLogger("pyturbo") logger.setLevel("DEBUG") stream_handler = StreamHandler() stream_handler.setLevel("DEBUG") - handler_format = Formatter('%(name)s - %(levelname)s - %(lineno)d - %(message)s') + handler_format = Formatter( + "%(name)s - %(levelname)s - %(lineno)d - %(message)s" + ) stream_handler.setFormatter(handler_format) logger.addHandler(stream_handler) - downloader_test_dir=os.path.join(pyturbo_root, "tests", "downloader") + downloader_test_dir = os.path.join(pyturbo_root, "tests", "downloader") os.chdir(downloader_test_dir) - # moved to examples \ No newline at end of file + # moved to examples diff --git a/turbogenius/pyturbo/utils/env.py b/turbogenius/pyturbo/utils/env.py index 196b31f..cdf6e36 100755 --- a/turbogenius/pyturbo/utils/env.py +++ b/turbogenius/pyturbo/utils/env.py @@ -13,21 +13,25 @@ from __future__ import print_function # python modules -import os, sys +import os import subprocess # set logger -from logging import config, getLogger, StreamHandler, Formatter -logger = getLogger('pyturbo').getChild(__name__) +from logging import getLogger -# pyturbo module -sys.path.append(os.path.dirname(os.path.abspath(__file__))) +logger = getLogger("pyturbo").getChild(__name__) # turbo-genius related path lists -pyturbo_root=os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(__file__)), '../../')) -pyturbo_source_dir=os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(__file__)), '../')) -pyturbo_data_dir=os.path.join(pyturbo_source_dir, 'data') -pyturbo_tmp_dir=os.path.join(os.path.abspath(os.environ['HOME']), '.pyturbo_tmp') +pyturbo_root = os.path.abspath( + os.path.join(os.path.dirname(os.path.abspath(__file__)), "../../") +) +pyturbo_source_dir = os.path.abspath( + os.path.join(os.path.dirname(os.path.abspath(__file__)), "../") +) +pyturbo_data_dir = os.path.join(pyturbo_source_dir, "data") +pyturbo_tmp_dir = os.path.join( + os.path.abspath(os.environ["HOME"]), ".pyturbo_tmp" +) # generate pyturbo temp. dir. os.makedirs(pyturbo_tmp_dir, exist_ok=True) @@ -40,20 +44,30 @@ try: sys_env = os.environ.copy() cmd = "which readalles.x" - turborvb_root = os.path.dirname(os.path.dirname(subprocess.check_output(cmd, shell=True, env=sys_env))).decode() + turborvb_root = os.path.dirname( + os.path.dirname( + subprocess.check_output(cmd, shell=True, env=sys_env) + ) + ).decode() turborvb_bin_root = os.path.join(turborvb_root, "bin") except: - raise ValueError("Set TURBORVB_ROOT (e.g., export TURBORVB_ROOT=XXX in ~.bashrc)") + raise ValueError( + "Set TURBORVB_ROOT (e.g., export TURBORVB_ROOT=XXX in ~.bashrc)" + ) if "TURBOMAKEFORT10_RUN_COMMAND" in os.environ: turbo_makefort10_run_command = os.environ["TURBOMAKEFORT10_RUN_COMMAND"] else: turbo_makefort10_run_command = "makefort10.x" if "TURBOCONVERTFORT10_RUN_COMMAND" in os.environ: - turbo_convertfort10_run_command = os.environ["TURBOCONVERTFORT10_RUN_COMMAND"] + turbo_convertfort10_run_command = os.environ[ + "TURBOCONVERTFORT10_RUN_COMMAND" + ] else: turbo_convertfort10_run_command = "convertfort10.x" if "TURBOCONVERTFORT10MOL_RUN_COMMAND" in os.environ: - turbo_convertfort10mol_run_command = os.environ["TURBOCONVERTFORT10MOL_RUN_COMMAND"] + turbo_convertfort10mol_run_command = os.environ[ + "TURBOCONVERTFORT10MOL_RUN_COMMAND" + ] else: turbo_convertfort10mol_run_command = "convertfort10mol.x" if "TURBOPREP_RUN_COMMAND" in os.environ: @@ -73,11 +87,15 @@ else: turbo_forcevmc_run_command = "forcevmc.sh" if "TURBOFORCEVMC_KPOINTS_RUN_COMMAND" in os.environ: - turbo_forcevmc_kpoints_run_command = os.environ["TURBOFORCEVMC_KPOINTS_RUN_COMMAND"] + turbo_forcevmc_kpoints_run_command = os.environ[ + "TURBOFORCEVMC_KPOINTS_RUN_COMMAND" + ] else: turbo_forcevmc_kpoints_run_command = "forcevmc_kpoints.sh" if "TURBOFORCEVMC_KPOINTS_PARA_RUN_COMMAND" in os.environ: - turbo_forcevmc_kpoints_para_run_command = os.environ["TURBOFORCEVMC_KPOINTS_PARA_RUN_COMMAND"] + turbo_forcevmc_kpoints_para_run_command = os.environ[ + "TURBOFORCEVMC_KPOINTS_PARA_RUN_COMMAND" + ] else: turbo_forcevmc_kpoints_para_run_command = "forcevmc_kpoints_parallel.sh" if "TURBOFORCEFN_RUN_COMMAND" in os.environ: @@ -85,11 +103,15 @@ else: turbo_forcefn_run_command = "forcefn.sh" if "TURBOFORCEFN_KPOINTS_RUN_COMMAND" in os.environ: - turbo_forcefn_kpoints_run_command = os.environ["TURBOFORCEFN_KPOINTS_RUN_COMMAND"] + turbo_forcefn_kpoints_run_command = os.environ[ + "TURBOFORCEFN_KPOINTS_RUN_COMMAND" + ] else: turbo_forcefn_kpoints_run_command = "forcefn_kpoints.sh" if "TURBOFORCEFN_KPOINTS_PARA_RUN_COMMAND" in os.environ: - turbo_forcefn_kpoints_para_run_command = os.environ["TURBOFORCEFN_KPOINTS_PARA_RUN_COMMAND"] + turbo_forcefn_kpoints_para_run_command = os.environ[ + "TURBOFORCEFN_KPOINTS_PARA_RUN_COMMAND" + ] else: turbo_forcefn_kpoints_para_run_command = "forcefn_kpoints_parallel.sh" if "TURBOCOPYJAS_RUN_COMMAND" in os.environ: @@ -97,6 +119,8 @@ else: turbo_copyjas_command = "copyjas.x" if "TURBOCONVERTPFAFF_RUN_COMMAND" in os.environ: - turbo_convertfortpfaff_run_command = os.environ["TURBOCONVERTPFAFF_RUN_COMMAND"] + turbo_convertfortpfaff_run_command = os.environ[ + "TURBOCONVERTPFAFF_RUN_COMMAND" + ] else: turbo_convertfortpfaff_run_command = "convertpfaff.x" diff --git a/turbogenius/pyturbo/utils/execute.py b/turbogenius/pyturbo/utils/execute.py index 4b49da6..7bae15d 100755 --- a/turbogenius/pyturbo/utils/execute.py +++ b/turbogenius/pyturbo/utils/execute.py @@ -13,43 +13,54 @@ from __future__ import print_function # python modules -import os, sys +import os +import sys import re import subprocess +from typing import Optional # set logger -from logging import config, getLogger, StreamHandler, Formatter -logger = getLogger('pyturbo').getChild(__name__) +from logging import getLogger -# turbogenius module -sys.path.append(os.path.dirname(os.path.abspath(__file__))) +logger = getLogger("pyturbo").getChild(__name__) -def run(binary, input_name = None, output_name = "out.o"): + +def run( + binary: str, input_name: Optional[str] = None, output_name: str = "out.o" +): sys_env = os.environ.copy() if input_name is None: cmd = f"{binary} > {output_name}" else: cmd = f"{binary} < {input_name} > {output_name}" - #p = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True, check=True, env=my_env) - #logger.info(p.stdout.decode()) - #logger.info(p.stderr.decode()) + # p = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True, check=True, env=my_env) + # logger.info(p.stdout.decode()) + # logger.info(p.stderr.decode()) logger.info(f"Execute command(s): {cmd}") - if sys.platform == 'darwin': + if sys.platform == "darwin": if "LD_LIBRARY_PATH" in sys_env: - if re.match(r".*bash.*", sys_env["SHELL"]) or re.match(r".*zsh.*", sys_env["SHELL"]): + if re.match(r".*bash.*", sys_env["SHELL"]) or re.match( + r".*zsh.*", sys_env["SHELL"] + ): cmd = f"export LD_LIBRARY_PATH={sys_env['LD_LIBRARY_PATH']} && {cmd}" - elif re.match(r".*csh.*", sys_env["SHELL"]) or re.match(r".*tcsh.*", sys_env["SHELL"]): + elif re.match(r".*csh.*", sys_env["SHELL"]) or re.match( + r".*tcsh.*", sys_env["SHELL"] + ): cmd = f"setenv LD_LIBRARY_PATH {sys_env['LD_LIBRARY_PATH']} && {cmd}" else: raise NotImplementedError if "DYLD_LIBRARY_PATH" in sys_env: - if re.match(r".*bash.*", sys_env["SHELL"]) or re.match(r".*zsh.*", sys_env["SHELL"]): + if re.match(r".*bash.*", sys_env["SHELL"]) or re.match( + r".*zsh.*", sys_env["SHELL"] + ): cmd = f"export DYLD_LIBRARY_PATH={sys_env['DYLD_LIBRARY_PATH']} && {cmd}" - elif re.match(r".*csh.*", sys_env["SHELL"]) or re.match(r".*tcsh.*", sys_env["SHELL"]): + elif re.match(r".*csh.*", sys_env["SHELL"]) or re.match( + r".*tcsh.*", sys_env["SHELL"] + ): cmd = f"setenv DYLD_LIBRARY_PATH {sys_env['DYLD_LIBRARY_PATH']} && {cmd}" else: raise NotImplementedError - subprocess.check_call(cmd, shell=True, env=sys_env) \ No newline at end of file + subprocess.check_call(cmd, shell=True, env=sys_env) diff --git a/turbogenius/pyturbo/utils/units.py b/turbogenius/pyturbo/utils/units.py index 644686d..ea59785 100755 --- a/turbogenius/pyturbo/utils/units.py +++ b/turbogenius/pyturbo/utils/units.py @@ -10,15 +10,12 @@ from __future__ import print_function # python modules -import os -import sys + # set logger -from logging import config, getLogger, StreamHandler, Formatter -logger = getLogger('pyturbo').getChild(__name__) +from logging import getLogger -# turbogenius module -sys.path.append(os.path.dirname(os.path.abspath(__file__))) +logger = getLogger("pyturbo").getChild(__name__) """ Units used in Turbo-Genius @@ -29,11 +26,10 @@ are defined as 1.0 """ -#Length +# Length Bohr = 1.0 -Angstrom = 1.0 / 0.529177210903 # Bohr +Angstrom = 1.0 / 0.529177210903 # Bohr -#Energy +# Energy Ha = 1.0 -Ry = 2.0 # Ha - +Ry = 2.0 # Ha diff --git a/turbogenius/pyturbo/utils/utility.py b/turbogenius/pyturbo/utils/utility.py index 93a21a6..b84e00d 100755 --- a/turbogenius/pyturbo/utils/utility.py +++ b/turbogenius/pyturbo/utils/utility.py @@ -96,8 +96,8 @@ def prompt(text, checker): def get_str_variable_type_auto(variable): - logger.debug(f"variable={variable}") - logger.debug(f"isdecimal={variable.isdecimal()}") + # logger.debug(f"variable={variable}") + # logger.debug(f"isdecimal={variable.isdecimal()}") if not variable.replace("-", "").isdecimal(): try: diff --git a/turbogenius/pyturbo/vmc.py b/turbogenius/pyturbo/vmc.py index 346a1f6..a098a92 100755 --- a/turbogenius/pyturbo/vmc.py +++ b/turbogenius/pyturbo/vmc.py @@ -17,6 +17,7 @@ import os import re import numpy as np +from typing import Optional # turbo-genius modules from turbogenius.pyturbo.namelist import Namelist @@ -46,10 +47,12 @@ class VMC(FortranIO): def __init__( self, - in_fort10="fort.10", - namelist=Namelist(), - twist_average=False, # False or 0: single-k, True or 1: Monkhorst-Pack, 2: manual k-grid + in_fort10: str = "fort.10", + namelist: Optional[Namelist] = None, + twist_average: bool = False, # False or 0: single-k, True or 1: Monkhorst-Pack, 2: manual k-grid ): + if namelist is None: + namelist = Namelist() """ input values @@ -99,7 +102,7 @@ def __str__(self): def sanity_check(self): pass - def generate_input(self, input_name="datasvmc.input"): + def generate_input(self, input_name: str = "datasvmc.input"): self.namelist.write(input_name) # check if twist_average is manual if self.twist_average == 2: # k-points are set manually @@ -137,7 +140,9 @@ def generate_input(self, input_name="datasvmc.input"): f.writelines(lines) logger.info(f"{input_name} has been generated.") - def run(self, input_name="datasvmc.input", output_name="out_vmc"): + def run( + self, input_name: str = "datasvmc.input", output_name: str = "out_vmc" + ): remove_file(file="pip0.d") remove_file(file="forces.dat") run( @@ -146,7 +151,9 @@ def run(self, input_name="datasvmc.input", output_name="out_vmc"): output_name=output_name, ) - def check_results(self, output_names=["out_vmc"]): + def check_results(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_vmc"] flags = [] for output_name in output_names: file_check(output_name) @@ -160,8 +167,11 @@ def check_results(self, output_names=["out_vmc"]): flags.append(False) return flags - def get_estimated_time_for_1_generation(self, output_names=["out_vmc"]): - + def get_estimated_time_for_1_generation( + self, output_names: Optional[list] = None + ): + if output_names is None: + output_names = ["out_vmc"] out_min = [] for output_name in output_names: with open(output_name, "r") as f: @@ -186,7 +196,7 @@ def get_estimated_time_for_1_generation(self, output_names=["out_vmc"]): return ave_time_1_generation # sec. @staticmethod - def read_energy(twist_average=False): + def read_energy(twist_average: bool = False): if twist_average: line = get_line_from_file(file="pip0.d", line_no=0).split() energy = float(line[3]) @@ -197,7 +207,13 @@ def read_energy(twist_average=False): error = float(line[3]) return energy, error - def get_energy(self, init=10, bin=10, num_proc=-1, rerun=False): + def get_energy( + self, + init: int = 10, + bin: int = 10, + num_proc: int = -1, + rerun: bool = False, + ): force_compute_flag = False if rerun: force_compute_flag = True @@ -224,7 +240,13 @@ def get_energy(self, init=10, bin=10, num_proc=-1, rerun=False): logger.debug("energy={}, error={}".format(energy, error)) return energy, error - def get_forces(self, init=10, bin=10, num_proc=-1, rerun=False): + def get_forces( + self, + init: int = 10, + bin: int = 10, + num_proc: int = -1, + rerun: bool = False, + ): force_compute_flag = False if rerun: @@ -324,7 +346,9 @@ def get_forces(self, init=10, bin=10, num_proc=-1, rerun=False): # unit (Ha/au) return force_matrix, force_matrix_error_bar - def compute_energy_and_forces(self, init=10, bin=10, pulay=1, num_proc=-1): + def compute_energy_and_forces( + self, init: int = 10, bin: int = 10, pulay: int = 1, num_proc: int = -1 + ): if self.twist_average: if num_proc == -1: logger.warning( @@ -336,7 +360,7 @@ def compute_energy_and_forces(self, init=10, bin=10, pulay=1, num_proc=-1): num_proc = os.cpu_count() if num_proc > 1: command = turbo_forcevmc_kpoints_para_run_command - #command = turbo_forcevmc_kpoints_run_command # for the time being!!! because paperoga does not have sufficient memory. + # command = turbo_forcevmc_kpoints_run_command # for the time being!!! because paperoga does not have sufficient memory. else: command = turbo_forcevmc_kpoints_run_command cmd = "{:s} {:d} {:d} {:d} {:d}".format( @@ -358,13 +382,13 @@ def read_default_namelist(): return namelist @staticmethod - def read_namelist_from_file(file): + def read_namelist_from_file(file: str): namelist = Namelist.parse_namelist_from_file(file) return namelist @classmethod def parse_from_default_namelist( - cls, in_fort10="fort.10", twist_average=False + cls, in_fort10: str = "fort.10", twist_average: bool = False ): namelist = cls.read_default_namelist() return cls( @@ -372,7 +396,9 @@ def parse_from_default_namelist( ) @classmethod - def parse_from_file(cls, file, in_fort10="fort.10", twist_average=False): + def parse_from_file( + cls, file: str, in_fort10: str = "fort.10", twist_average: bool = False + ): namelist = Namelist.parse_namelist_from_file(file) return cls( in_fort10=in_fort10, namelist=namelist, twist_average=twist_average diff --git a/turbogenius/pyturbo/vmcopt.py b/turbogenius/pyturbo/vmcopt.py index 0af6c31..133f410 100755 --- a/turbogenius/pyturbo/vmcopt.py +++ b/turbogenius/pyturbo/vmcopt.py @@ -21,6 +21,7 @@ import pandas as pd import matplotlib.pyplot as plt import matplotlib.ticker as ticker +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -50,10 +51,12 @@ class VMCopt(FortranIO): def __init__( self, - in_fort10="fort.10", - namelist=Namelist(), - twist_average=False, # False or 0: single-k, True or 1: Monkhorst-Pack, 2: manual k-grid + in_fort10: str = "fort.10", + namelist: Optional[Namelist] = None, + twist_average: bool = False, # False or 0: single-k, True or 1: Monkhorst-Pack, 2: manual k-grid ): + if namelist is None: + namelist = Namelist() """ input values @@ -86,12 +89,10 @@ def manual_kpoints(self, kpoints): parameter="yes_kpoints", value=".true.", namelist="¶meters" ) # self.namelist.set_parameter(parameter="yeswrite10", value=".true.", namelist="&optimization") - self.namelist.set_parameter( - parameter="kp_type", value=2, namelist="&kpoints" - ) + self.namelist.set_parameter(parameter="kp_type", value=2, namelist="&kpoints") self.namelist.set_parameter( parameter="nk1", - value=len(kpoints_up) + len(kpoints_dn), + value=len(kpoints_up), namelist="&kpoints", ) @@ -105,7 +106,7 @@ def __str__(self): def sanity_check(self): pass - def generate_input(self, input_name): + def generate_input(self, input_name: str = "datasmin.input"): self.namelist.write(input_name) # check if twist_average is manual if self.twist_average == 2: # k-points are set manually @@ -114,8 +115,7 @@ def generate_input(self, input_name): with open(input_name, "r") as f: lines = f.readlines() kpoint_index = [ - True if re.match(r".*&kpoints.*", line) else False - for line in lines + True if re.match(r".*&kpoints.*", line) else False for line in lines ].index(True) insert_lineno = -1 for i, line in enumerate(lines[kpoint_index + 1 :]): @@ -143,7 +143,7 @@ def generate_input(self, input_name): f.writelines(lines) logger.info(f"{os.path.basename(input_name)} has been generated.") - def run(self, input_name="datasmin.input", output_name="out_min"): + def run(self, input_name: str = "datasmin.input", output_name: str = "out_min"): run( turbo_qmc_run_command, input_name=input_name, @@ -151,24 +151,25 @@ def run(self, input_name="datasmin.input", output_name="out_min"): ) remove_file(file="stop.dat") - def check_results(self, output_names=["out_min"]): + def check_results(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_min"] flags = [] for output_name in output_names: file_check(output_name) with open(output_name, "r") as f: lines = f.readlines() - if any( - [re.match(r".*Final.*tstep.*found.*", line) for line in lines] - ): + if any([re.match(r".*Final.*tstep.*found.*", line) for line in lines]): flags.append(True) else: flags.append(False) return flags def plot_energy_and_devmax( - self, output_names=["out_min"], interactive=False + self, output_names: Optional[list] = None, interactive: bool = False ): - + if output_names is None: + output_names = ["out_min"] # plot the energies and devmax logger.info("Plotting the energies and devmax") out_min = [] @@ -189,11 +190,7 @@ def plot_energy_and_devmax( devmax_list = list( map( lambda x: x.split(), - [ - i - for i in out_min - if re.match(r".*devmax.*par.*Normal.*", i) - ], + [i for i in out_min if re.match(r".*devmax.*par.*Normal.*", i)], ) ) devmax_pandas_str = pd.DataFrame(devmax_list, columns=col) @@ -251,13 +248,12 @@ def plot_energy_and_devmax( ax1.legend(h1 + h2, l1 + l2, loc="upper right") if interactive: plt.waitforbuttonpress() - plt.savefig( - "plot_energy_and_devmax.png", bbox_inches="tight", pad_inches=0.2 - ) + plt.savefig("plot_energy_and_devmax.png", bbox_inches="tight", pad_inches=0.2) plt.close() - def get_energy(self, output_names=["out_min"]): - + def get_energy(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_min"] out_min = [] for output_name in output_names: with open(output_name, "r") as f: @@ -266,28 +262,21 @@ def get_energy(self, output_names=["out_min"]): energy_list = list( map( float, - [ - i.split()[3] - for i in out_min - if re.match(r".*New.*Energy.*", i) - ], + [i.split()[3] for i in out_min if re.match(r".*New.*Energy.*", i)], ) ) error_list = list( map( float, - [ - i.split()[4] - for i in out_min - if re.match(r".*New.*Energy.*", i) - ], + [i.split()[4] for i in out_min if re.match(r".*New.*Energy.*", i)], ) ) return energy_list, error_list - def get_devmax(self, output_names=["out_min"]): - + def get_devmax(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_min"] out_min = [] for output_name in output_names: with open(output_name, "r") as f: @@ -306,7 +295,9 @@ def get_devmax(self, output_names=["out_min"]): return devmax_list - def get_estimated_time_for_1_generation(self, output_names=["out_min"]): + def get_estimated_time_for_1_generation(self, output_names: Optional[list] = None): + if output_names is None: + output_names = ["out_min"] out_min = [] for output_name in output_names: @@ -319,19 +310,15 @@ def get_estimated_time_for_1_generation(self, output_names=["out_min"]): [ i.split()[5] for i in out_min - if re.match( - r".*Average.*time.*for.*1000.*generations.*", i - ) + if re.match(r".*Average.*time.*for.*1000.*generations.*", i) ], ) ) - ave_time_1_generation = ( - np.mean(np.array(ave_time_1000_generations)) / 1000 - ) + ave_time_1_generation = np.mean(np.array(ave_time_1000_generations)) / 1000 return ave_time_1_generation # sec. - def plot_parameters_history(self, interactive=True): + def plot_parameters_history(self, interactive: bool = True): # save parameters current_dir = os.getcwd() cmd = f"(echo '1 1 0 0'; echo '0'; echo '100000') | {os.path.join(turborvb_bin_root, 'readalles.x')}" @@ -387,9 +374,7 @@ def plot_parameters_history(self, interactive=True): logger.info("KeyboardInterrupt") break plt.savefig( - os.path.join( - graph_save_dir, "Parameter_No{}_all.png".format(i) - ), + os.path.join(graph_save_dir, "Parameter_No{}_all.png".format(i)), bbox_inches="tight", pad_inches=0.2, ) @@ -397,9 +382,9 @@ def plot_parameters_history(self, interactive=True): def average_optimized_parameters( self, - equil_steps=10, - input_file_used="datasmin.input", - graph_plot=False, + equil_steps: int = 10, + input_file_used: str = "datasmin.input", + graph_plot: bool = False, ): """ @@ -469,9 +454,7 @@ def average_optimized_parameters( linestyle="dashed", label="averaged", ) - plt.xlabel( - "Iteration", fontname="Times New Roman", fontsize=14 - ) + plt.xlabel("Iteration", fontname="Times New Roman", fontsize=14) plt.ylabel("Value", fontname="Times New Roman", fontsize=14) # plt.legend(frameon=True) plt.title("Parameter_No.{}".format(i)) @@ -497,11 +480,11 @@ def average_optimized_parameters( vmcopt = VMCopt.parse_from_file( file=input_file_used, in_fort10=self.in_fort10, - twist_average=self.twist_average, + twist_average=False, # always False because we do not need twist for the average. ) vmcopt.set_parameter("iopt", 1) vmcopt.set_parameter("ngen", 0) - if self.twist_average: + if self.twist_average != 0: vmcopt.set_parameter("yes_kpoints", ".false.", "¶meters") # replace iopt from 0 -> 1 @@ -544,32 +527,28 @@ def average_optimized_parameters( @staticmethod def read_default_namelist(): - vmcopt_default_file = os.path.join( - pyturbo_data_dir, "vmcopt", "datasmin.input" - ) + vmcopt_default_file = os.path.join(pyturbo_data_dir, "vmcopt", "datasmin.input") namelist = Namelist.parse_namelist_from_file(vmcopt_default_file) return namelist @staticmethod - def read_namelist_from_file(file): + def read_namelist_from_file(file: str): namelist = Namelist.parse_namelist_from_file(file) return namelist @classmethod def parse_from_default_namelist( - cls, in_fort10="fort.10", twist_average=False + cls, in_fort10: str = "fort.10", twist_average: bool = False ): namelist = cls.read_default_namelist() - return cls( - in_fort10=in_fort10, namelist=namelist, twist_average=twist_average - ) + return cls(in_fort10=in_fort10, namelist=namelist, twist_average=twist_average) @classmethod - def parse_from_file(cls, file, in_fort10="fort.10", twist_average=False): + def parse_from_file( + cls, file, in_fort10: str = "fort.10", twist_average: bool = False + ): namelist = Namelist.parse_namelist_from_file(file) - return cls( - in_fort10=in_fort10, namelist=namelist, twist_average=twist_average - ) + return cls(in_fort10=in_fort10, namelist=namelist, twist_average=twist_average) if __name__ == "__main__": @@ -577,9 +556,7 @@ def parse_from_file(cls, file, in_fort10="fort.10", twist_average=False): logger.setLevel("INFO") stream_handler = StreamHandler() stream_handler.setLevel("DEBUG") - handler_format = Formatter( - "%(name)s - %(levelname)s - %(lineno)d - %(message)s" - ) + handler_format = Formatter("%(name)s - %(levelname)s - %(lineno)d - %(message)s") stream_handler.setFormatter(handler_format) logger.addHandler(stream_handler) diff --git a/turbogenius/readforward_genius.py b/turbogenius/readforward_genius.py index 694a58f..c9ba3f1 100755 --- a/turbogenius/readforward_genius.py +++ b/turbogenius/readforward_genius.py @@ -12,6 +12,7 @@ # python modules import os +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -114,7 +115,7 @@ def run( flags = self.readforward.check_results(output_names=[output_name]) assert all(flags) - def check_results(self, output_names: list = ["out_readforward"]) -> None: + def check_results(self, output_names: Optional[list] = None) -> bool: """ Check the result. @@ -123,6 +124,8 @@ def check_results(self, output_names: list = ["out_readforward"]) -> None: Returns: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_readforward"] return self.readforward.check_results(output_names=output_names) diff --git a/turbogenius/tools_genius.py b/turbogenius/tools_genius.py index 8d6d43f..13feb8c 100755 --- a/turbogenius/tools_genius.py +++ b/turbogenius/tools_genius.py @@ -30,7 +30,9 @@ def copy_jastrow( - fort10_to="fort.10", fort10_from="fort.10_new", twist_flag=False + fort10_to: str = "fort.10", + fort10_from: str = "fort.10_new", + twist_flag: bool = False, ): """ Copy Jastrow factors diff --git a/turbogenius/trexio_to_turborvb.py b/turbogenius/trexio_to_turborvb.py index 4ac8bac..65d0b83 100644 --- a/turbogenius/trexio_to_turborvb.py +++ b/turbogenius/trexio_to_turborvb.py @@ -19,10 +19,14 @@ import argparse import numpy as np import glob +from typing import Optional # logger from logging import getLogger, StreamHandler, Formatter +# import trexio +import trexio + # import pyturbo modules from turbogenius.pyturbo.structure import Structure, Cell from turbogenius.pyturbo.pseudopotentials import Pseudopotentials @@ -53,7 +57,7 @@ def trexio_to_turborvb_wf( trexio_file: str, - jas_basis_sets=Jas_Basis_sets(), + jas_basis_sets: Optional[Jas_Basis_sets] = None, max_occ_conv: int = 0, mo_num_conv: int = -1, only_mol: bool = True, @@ -70,6 +74,8 @@ def trexio_to_turborvb_wf( only_mol (bool): if True, only moleculer orbitals option = True in convertfort10mol cleanup (bool): clean up temporary files """ + if jas_basis_sets is None: + jas_basis_sets = Jas_Basis_sets() # os.environ["DYLD_LIBRARY_PATH"]=os.environ["LIBRARY_PATH"] # prefix and file names @@ -139,7 +145,6 @@ def trexio_to_turborvb_wf( phase_dn = [0.0, 0.0, 0.0] complex_flag = trexio_r.complex_flag - logger.info(complex_flag) # view # from ase.visualize import view # atom=structure.get_ase_atom() @@ -164,14 +169,18 @@ def trexio_to_turborvb_wf( logger.debug(basis_shell_index) # Pseudo potentials info - ecp_max_ang_mom_plus_1 = trexio_r.ecp_max_ang_mom_plus_1 - ecp_z_core = trexio_r.ecp_z_core - ecp_num = trexio_r.ecp_num - ecp_ang_mom = trexio_r.ecp_ang_mom - ecp_nucleus_index = trexio_r.ecp_nucleus_index - ecp_exponent = trexio_r.ecp_exponent - ecp_coefficient = trexio_r.ecp_coefficient - ecp_power = trexio_r.ecp_power + try: + ecp_max_ang_mom_plus_1 = trexio_r.ecp_max_ang_mom_plus_1 + ecp_z_core = trexio_r.ecp_z_core + ecp_num = trexio_r.ecp_num + ecp_ang_mom = trexio_r.ecp_ang_mom + ecp_nucleus_index = trexio_r.ecp_nucleus_index + ecp_exponent = trexio_r.ecp_exponent + ecp_coefficient = trexio_r.ecp_coefficient + ecp_power = trexio_r.ecp_power + has_ecp = True + except AttributeError: + has_ecp = False # ao info ao_cartesian = trexio_r.ao_cartesian @@ -187,21 +196,93 @@ def trexio_to_turborvb_wf( mo_num = trexio_r.mo_num mo_coefficient = trexio_r.mo_coefficient mo_occupation = trexio_r.mo_occupation + mo_spin = trexio_r.mo_spin if complex_flag: mo_coefficient_imag = trexio_r.mo_coefficient_imag mo_coefficient = [ [complex(i, j) for i, j in zip(mo_real, mo_imag)] for mo_real, mo_imag in zip(mo_coefficient, mo_coefficient_imag) ] + if all([spin == 0 for spin in mo_spin]) or all( + [spin == 1 for spin in mo_spin] + ): + logger.info("MOs are spin-restricted (i.e., alpha==beta).") + spin_restricted = True + elif all([spin == 0 or spin == 1 for spin in mo_spin]): + logger.info("MOs are spin-unrestricted (i.e., alpha!=beta).") + spin_restricted = False + else: + raise ValueError + + # spin unrestricted is supported only with trexio >= 1.3.0 + if not spin_restricted and not trexio.__version__ >= '1.3.0': + logger.error("spin unrestricted is supported only with trexio >= 1.3.0") + raise NotImplementedError - logger.debug(ecp_max_ang_mom_plus_1) - logger.debug(ecp_z_core) - logger.debug(ecp_num) - logger.debug(ecp_ang_mom) - logger.debug(ecp_nucleus_index) - logger.debug(ecp_exponent) - logger.debug(ecp_coefficient) - logger.debug(ecp_power) + # check if the num. of MOs for alpha and beta are the same. + if not spin_restricted: + if list(mo_spin).count(0) != list(mo_spin).count(1): + logger.error( + "The number of alpha- and beta-MOs are not consistent!!" + ) + raise ValueError + + # At present, this python script assumes that MOs are [alpha, alpha..... alpha, beta....beta] for + # an unrestriced WF. Check it. + if not spin_restricted: + if not ( + all([spin == 0 for spin in mo_spin[0 : int(len(mo_spin) / 2)]]) + and all( + [ + spin == 1 + for spin in mo_spin[int(len(mo_spin) / 2) : len(mo_spin)] + ] + ) + ): + logger.error( + "The MOs ordering is not [alpha, ... alpha, beta, .... beta]" + ) + raise NotImplementedError + + # At present, this python script assumes that the occupation is in the descending order. check it. + if spin_restricted: + if not all( + [ + mo_occ == mo_occ_sorted + for mo_occ, mo_occ_sorted in zip( + mo_occupation, sorted(mo_occupation, reverse=True) + ) + ] + ): + raise NotImplementedError + else: + # alpha spin + if not all( + [ + mo_occ == mo_occ_sorted + for mo_occ, mo_occ_sorted in zip( + mo_occupation[0 : int(len(mo_spin) / 2)], + sorted( + mo_occupation[0 : int(len(mo_spin) / 2)], reverse=True + ), + ) + ] + ): + raise NotImplementedError + # beta spin + if not all( + [ + mo_occ == mo_occ_sorted + for mo_occ, mo_occ_sorted in zip( + mo_occupation[int(len(mo_spin) / 2) : len(mo_spin)], + sorted( + mo_occupation[int(len(mo_spin) / 2) : len(mo_spin)], + reverse=True, + ), + ) + ] + ): + raise NotImplementedError # basis sets shell_ang_mom_turbo_notation = [] @@ -234,21 +315,24 @@ def trexio_to_turborvb_wf( ) # Pseudopotentials - cutoff = [0.0] * len(ecp_z_core) - pseudopotentials = Pseudopotentials( - max_ang_mom_plus_1=ecp_max_ang_mom_plus_1, - z_core=ecp_z_core, - cutoff=cutoff, - nucleus_index=ecp_nucleus_index, - element_list=element_list, - ang_mom=ecp_ang_mom, - exponent=ecp_exponent, - coefficient=ecp_coefficient, - power=ecp_power, - ) - pseudopotentials.set_cutoffs() - logger.debug(pseudopotentials.cutoff) - pseudopotentials.write_pseudopotential_turborvb_file() + if has_ecp: + cutoff = [0.0] * len(ecp_z_core) + pseudopotentials = Pseudopotentials( + max_ang_mom_plus_1=ecp_max_ang_mom_plus_1, + z_core=ecp_z_core, + cutoff=cutoff, + nucleus_index=ecp_nucleus_index, + element_list=element_list, + ang_mom=ecp_ang_mom, + exponent=ecp_exponent, + coefficient=ecp_coefficient, + power=ecp_power, + ) + pseudopotentials.set_cutoffs() + logger.debug(pseudopotentials.cutoff) + pseudopotentials.write_pseudopotential_turborvb_file() + else: + pseudopotentials = Pseudopotentials() # set jastrow_type: if jas_basis_sets.shell_num == 0: @@ -316,11 +400,23 @@ def trexio_to_turborvb_wf( parameter="phasedo(3)", value=phase_dn[2], namelist="&system" ) + # spin-restricted or spin-unrestricted + # => + # symmetrized AGP (AGPs) or unsymmetrized AGP (AGPu) + if spin_restricted: + namelist.set_parameter( + parameter="symmagp", value=".true.", namelist="&symmetries" + ) + else: + namelist.set_parameter( + parameter="symmagp", value=".false.", namelist="&symmetries" + ) + + # complex or real if complex_flag: namelist.set_parameter( parameter="complexfort10", value=".true.", namelist="&system" ) - # logger.error("complex is not implemented in the converter.") else: namelist.set_parameter( parameter="complexfort10", value=".false.", namelist="&system" @@ -333,47 +429,115 @@ def trexio_to_turborvb_wf( pseudo_potentials=pseudopotentials, namelist=namelist, ) - makefort10.generate_input("makefort10.input") + makefort10.generate_input( + input_name="makefort10.input", basis_sets_unique_element=False + ) makefort10.run() # rename fort.10 shutil.move("fort.10_new", "fort.10_in") # here, specify how many MOs are used for the conversion!! - if max_occ_conv == 0 and mo_num_conv == -1: # default - logger.info(f"All the MOs ={mo_num} are used for the conversion.") - mo_num_use = mo_num - mo_coefficient = mo_coefficient - mo_occ = mo_occupation - elif max_occ_conv != 0 and mo_num_conv != -1: + + # note that the meaning of mo_num is different between ROHF and UHF calculations. + # ROHF -> MOs should have the same coeffs. for up and dn. This means that + # the unpaired MOs are included ONLY for the up part. (as it should be listed in the + # end of the MO list. In other words, the dn-pairs of the unpaired MOs (up) are not included in the + # beta part. On the other hand, in UHF cases, the dn-pairs of the unpaired MOs (up) + # can be included also in the beta part. + + # Anyway, here "mo_num_use" is the number of occupied, empty, and unpaired MOs for "alpha spin". + + if max_occ_conv != 0 and mo_num_conv != -1: logger.error( "max_occ_conv and mo_num_conv options cannot be used at the same time." ) raise ValueError + + elif max_occ_conv == 0 and mo_num_conv == -1: # default + if spin_restricted: + logger.info(f"All the MOs = {mo_num} are used for the conversion.") + mo_num_use = mo_num + mo_coefficient = mo_coefficient + mo_spin = mo_spin + else: + logger.info(f"All the MOs = {mo_num} are used for the conversion.") + logger.warning("This is an spin-unrestricted WF.") + logger.warning( + f"Therefore, the num. of MOs used for alpha (beta) is {int(mo_num / 2)}" + ) + mo_num_use = int(mo_num / 2) + mo_coefficient = mo_coefficient + mo_spin = mo_spin + else: if max_occ_conv != 0: mo_used_index = [] - logger.info(f"mo_occ_thr={max_occ_conv}") - # logger.info(f"mo_occupation") - # logger.info(mo_occupation) - for i in range(len(mo_occupation)): - mo_used_index.append(i) - if i == len(mo_occupation) - 1: - break - if ( - mo_occupation[i] >= max_occ_conv - and mo_occupation[i + 1] < max_occ_conv - ): - logger.info(f"mo_occ < {max_occ_conv} is 1-{i + 1}") - break - mo_num_use = len(mo_used_index) + if spin_restricted: + logger.info(f"mo_occ_thr={max_occ_conv}.") + for i in range(len(mo_occupation)): + mo_used_index.append(i) + if i == len(mo_occupation) - 1: + break + if ( + mo_occupation[i] >= max_occ_conv + and mo_occupation[i + 1] < max_occ_conv + ): + logger.info(f"mo_occ < {max_occ_conv} is 1-{i + 1}") + break + mo_num_use = len(mo_used_index) + else: # spin-unresticted + logger.info( + f"mo_occ_thr={max_occ_conv} both for alpha and beta spins." + ) + # spin-up (alpha) + for i in range(0, int(len(mo_occupation) / 2)): + mo_used_index.append(i) # alpha spin + if i == int(len(mo_occupation) / 2) - 1: + break + if ( + mo_occupation[i] >= max_occ_conv + and mo_occupation[i + 1] < max_occ_conv + ): + logger.info( + f"mo_occ(alpha) < {max_occ_conv} is 1-{i + 1}" + ) + break + # spin-dn (beta) + mo_used_index = mo_used_index + [ + i + int(len(mo_occupation) / 2) for i in mo_used_index + ] + + mo_num_use = int(len(mo_used_index) / 2) + mo_coefficient = [mo_coefficient[ppp] for ppp in mo_used_index] - mo_occ = [mo_occupation[ppp] for ppp in mo_used_index] + mo_occupation = [mo_occupation[ppp] for ppp in mo_used_index] + mo_spin = [mo_spin[ppp] for ppp in mo_used_index] else: # mo_num_conv != -1: - mo_num_use = mo_num_conv - mo_coefficient = mo_coefficient[0:mo_num_use] - mo_occ = mo_occupation[0:mo_num_use] + if spin_restricted: + mo_num_use = mo_num_conv + mo_coefficient = mo_coefficient[0:mo_num_use] + mo_occupation = mo_occupation[0:mo_num_use] + mo_spin = mo_spin[0:mo_num_use] + else: # spin-unrestricted + mo_num_use = mo_num_conv + mo_coefficient = ( + mo_coefficient[0:mo_num_use] + + mo_coefficient[ + int(mo_num / 2) : int(mo_num / 2) + mo_num_use + ] + ) + mo_occupation = ( + mo_occupation[0:mo_num_use] + + mo_occupation[ + int(mo_num / 2) : int(mo_num / 2) + mo_num_use + ] + ) + mo_spin = ( + mo_spin[0:mo_num_use] + + mo_spin[int(mo_num / 2) : int(mo_num / 2) + mo_num_use] + ) # in_convertfort10mol class convertfort10mol = Convertfort10mol.parse_from_default_namelist() @@ -396,492 +560,557 @@ def trexio_to_turborvb_wf( # check cartesian / spherical flag: if ao_cartesian != 0: - logger.debug("basis set is represent with cartesian") + logger.error("basis set is represent with cartesian") + raise NotImplementedError else: - logger.debug("basis set is represent with spherical") + logger.info("basis set is represent with spherical") mo_coefficient_turbo = [] mo_exponent_turbo = [] - num = mo_num_use - # num=1 - for mo_i in range(num): - logger.debug(f"============{mo_i+1}-th MO================") - mo_coefficient_list = [] - mo_exponent_list = [] - mo_nucleus_index_list = [] - mo_l_list = [] - mo_m_list = [] - # logger.debug(f"i={i}") - # logger.debug(mo_coefficient[i]) - # logger.debug(ao_shell) - - for ao_i, shell_index in enumerate(ao_shell): - logger.debug(f"--ao_i={ao_i}, shell_index={shell_index}--") - # for real spherical harmonic notations, - # see [https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics] - - # initialization - if ao_i == 0: - mo_coefficient_list_for_reordering = [] - mo_exponent_list_for_reordering = [] - mo_nucleus_index_list_for_reordering = [] - current_ang_mom = -1 - - # read ang_mom (i.e., angular momentum of the shell) - ang_mom = basis_shell_ang_mom[shell_index] - previous_ang_mom = current_ang_mom - current_ang_mom = ang_mom - - # set multiplicity - multiplicity = 2 * ang_mom + 1 - logger.debug(f"multiplicity = {multiplicity}") - - # check if the buffer is null, when ang_mom changes - if previous_ang_mom != current_ang_mom: - assert len(mo_coefficient_list_for_reordering) == 0 - assert len(mo_exponent_list_for_reordering) == 0 - - if current_ang_mom == 0: # s shell - logger.debug("s shell/no permutation is needed.") - logger.debug("(trexio notation): s(m=0)") - logger.debug("(makefun notation): s(m=0)") - reorder_index = [0] - reorder_m_list = [0] - reorder_l_list = [0] * 1 - - elif current_ang_mom == 1: # p shell - - if ao_cartesian != 0: - logger.debug("p shell/no permutation is needed.") - logger.debug( - "(trexio notation): px(m=+1), py(m=-1), pz(m=0)" - ) - logger.debug( - "(makefun notation): px(m=+1), py(m=-1), pz(m=0)" - ) - reorder_index = [0, 1, 2] - reorder_m_list = [+1, -1, 0] - reorder_l_list = [1] * 3 - else: - logger.debug("p shell/permutation is needed.") - logger.debug( - "(trexio notation): pz(m=0), px(m=+1), py(m=-1)" - ) - logger.debug( - "(makefun notation): px(m=+1), py(m=-1), pz(m=0)" - ) - reorder_index = [1, 2, 0] - reorder_m_list = [+1, -1, 0] - reorder_l_list = [1] * 3 - elif current_ang_mom == 2: # d shell + # Here is the most important part of the trexio_to_turborvb conversion. + # 1. Reordering the MOs. + # 2. Removing the duplicated exponents. because turbo does not compute them. - if ao_cartesian != 0: - logger.debug( - "Cartesian notation for d shell is not implemented yet! Sorry." - ) - raise NotImplementedError - else: - logger.debug("d shell/permutation is needed.") - logger.debug( - "(trexio notation): dz2(m=0), dzx=(m=+1), dyz(m=-1), dx2-y2(m=+2), dxy(m=-2)" - ) - logger.debug( - "(makefun notation): dz2(m=0), dx2-y2(m=+2), dxy(m=-2), dyz(m=-1), dzx=(m=+1)" - ) - reorder_index = [0, 3, 4, 2, 1] - reorder_m_list = [0, +2, -2, -1, +1] - reorder_l_list = [2] * 5 + # here one should remember that num := mo_num_use is the num for alpha (or beta) spin. + # for spin-restricted conversion, we should repeat this procedure twice. + if spin_restricted: + mo_i_shift_list = [0] + else: + mo_i_shift_list = [0, mo_num_use] + + for shift_v in mo_i_shift_list: + for mo_i in range(mo_num_use): + mo_i_shifted = mo_i + shift_v + logger.debug(f"============{mo_i_shifted+1}-th MO================") + mo_coefficient_list = [] + mo_exponent_list = [] + mo_nucleus_index_list = [] + mo_l_list = [] + mo_m_list = [] + # logger.debug(f"i={i}") + # logger.debug(mo_coefficient[i]) + # logger.debug(ao_shell) + + for ao_i, shell_index in enumerate(ao_shell): + logger.debug(f"--ao_i={ao_i}, shell_index={shell_index}--") + # for real spherical harmonic notations, + # see [https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics] + + # initialization + if ao_i == 0: + mo_coefficient_list_for_reordering = [] + mo_exponent_list_for_reordering = [] + mo_nucleus_index_list_for_reordering = [] + current_ang_mom = -1 + + # read ang_mom (i.e., angular momentum of the shell) + ang_mom = basis_shell_ang_mom[shell_index] + previous_ang_mom = current_ang_mom + current_ang_mom = ang_mom + + # set multiplicity + multiplicity = 2 * ang_mom + 1 + logger.debug(f"multiplicity = {multiplicity}") + + # check if the buffer is null, when ang_mom changes + if previous_ang_mom != current_ang_mom: + assert len(mo_coefficient_list_for_reordering) == 0 + assert len(mo_exponent_list_for_reordering) == 0 + + if current_ang_mom == 0: # s shell + logger.debug("s shell/no permutation is needed.") + logger.debug("(trexio notation): s(m=0)") + logger.debug("(makefun notation): s(m=0)") + reorder_index = [0] + reorder_m_list = [0] + reorder_l_list = [0] * 1 + + elif current_ang_mom == 1: # p shell + + if ao_cartesian != 0: + logger.debug("p shell/no permutation is needed.") + logger.debug( + "(trexio notation): px(m=+1), py(m=-1), pz(m=0)" + ) + logger.debug( + "(makefun notation): px(m=+1), py(m=-1), pz(m=0)" + ) + reorder_index = [0, 1, 2] + reorder_m_list = [+1, -1, 0] + reorder_l_list = [1] * 3 + else: + logger.debug("p shell/permutation is needed.") + logger.debug( + "(trexio notation): pz(m=0), px(m=+1), py(m=-1)" + ) + logger.debug( + "(makefun notation): px(m=+1), py(m=-1), pz(m=0)" + ) + reorder_index = [1, 2, 0] + reorder_m_list = [+1, -1, 0] + reorder_l_list = [1] * 3 - elif current_ang_mom == 3: # f shell + elif current_ang_mom == 2: # d shell - if ao_cartesian != 0: - logger.debug( - "Cartesian notation for f shell is not implemented yet! Sorry." - ) - raise NotImplementedError + if ao_cartesian != 0: + logger.debug( + "Cartesian notation for d shell is not implemented yet! Sorry." + ) + raise NotImplementedError + else: + logger.debug("d shell/permutation is needed.") + logger.debug( + "(trexio notation): dz2(m=0), dzx=(m=+1), dyz(m=-1), dx2-y2(m=+2), dxy(m=-2)" + ) + logger.debug( + "(makefun notation): dz2(m=0), dx2-y2(m=+2), dxy(m=-2), dyz(m=-1), dzx=(m=+1)" + ) + reorder_index = [0, 3, 4, 2, 1] + reorder_m_list = [0, +2, -2, -1, +1] + reorder_l_list = [2] * 5 - else: - logger.debug("f shell/no permutation is needed.") - logger.debug( - "(trexio notation): f3,0(l=0), f3,+1(l=+1), f3,-1(l=-1), f3,+2(l=+2), f3,-2(l=-2), f3,+3(l=+3), f3,-3(l=-3)" - ) - logger.debug( - "(makefun notation): f3,0(l=0), f3,+1(l=+1), f3,-1(l=-1), f3,+2(l=+2), f3,-2(l=-2), f3,+3(l=+3), f3,-3(l=-3)" - ) - reorder_index = [0, 1, 2, 3, 4, 5, 6] - reorder_m_list = [0, +1, -1, +2, -2, +3, -3] - reorder_l_list = [3] * 7 + elif current_ang_mom == 3: # f shell - elif current_ang_mom == 4: # g shell + if ao_cartesian != 0: + logger.debug( + "Cartesian notation for f shell is not implemented yet! Sorry." + ) + raise NotImplementedError - if ao_cartesian != 0: - logger.debug( - "Cartesian notation for g shell is not implemented yet! Sorry." - ) - raise NotImplementedError + else: + logger.debug("f shell/no permutation is needed.") + logger.debug( + "(trexio notation): f3,0(l=0), f3,+1(l=+1), f3,-1(l=-1), f3,+2(l=+2), f3,-2(l=-2), f3,+3(l=+3), f3,-3(l=-3)" + ) + logger.debug( + "(makefun notation): f3,0(l=0), f3,+1(l=+1), f3,-1(l=-1), f3,+2(l=+2), f3,-2(l=-2), f3,+3(l=+3), f3,-3(l=-3)" + ) + reorder_index = [0, 1, 2, 3, 4, 5, 6] + reorder_m_list = [0, +1, -1, +2, -2, +3, -3] + reorder_l_list = [3] * 7 - else: - logger.debug("g shell/no permutation is needed.") - logger.debug( - "(trexio notation): g4,0(l=0), g4,+1(l=+1), g4,-1(l=-1), g4,+2(l=+2), g4,-2(l=-2), g4,+3(l=+3), g4,-3(l=-3), g4,+4(l=+4), g4,-4(l=-4)" - ) - logger.debug( - "(makefun notation): g4,0(l=0), g4,+1(l=+1), g4,-1(l=-1), g4,+2(l=+2), g4,-2(l=-2), g4,+3(l=+3), g4,-3(l=-3), g4,+4(l=+4), g4,-4(l=-4)" - ) - reorder_index = [0, 1, 2, 3, 4, 5, 6, 7, 8] - reorder_m_list = [0, +1, -1, +2, -2, +3, -3, +4, -4] - reorder_l_list = [4] * 9 + elif current_ang_mom == 4: # g shell - elif current_ang_mom == 5: # h shell + if ao_cartesian != 0: + logger.debug( + "Cartesian notation for g shell is not implemented yet! Sorry." + ) + raise NotImplementedError - if ao_cartesian != 0: - logger.debug( - "Cartesian notation for h shell is not implemented yet! Sorry." - ) - raise NotImplementedError + else: + logger.debug("g shell/no permutation is needed.") + logger.debug( + "(trexio notation): g4,0(l=0), g4,+1(l=+1), g4,-1(l=-1), g4,+2(l=+2), g4,-2(l=-2), g4,+3(l=+3), g4,-3(l=-3), g4,+4(l=+4), g4,-4(l=-4)" + ) + logger.debug( + "(makefun notation): g4,0(l=0), g4,+1(l=+1), g4,-1(l=-1), g4,+2(l=+2), g4,-2(l=-2), g4,+3(l=+3), g4,-3(l=-3), g4,+4(l=+4), g4,-4(l=-4)" + ) + reorder_index = [0, 1, 2, 3, 4, 5, 6, 7, 8] + reorder_m_list = [0, +1, -1, +2, -2, +3, -3, +4, -4] + reorder_l_list = [4] * 9 - else: - logger.debug("h shell/no permutation is needed.") - logger.debug( - "(trexio notation): h5,0(l=0), h5,+1(l=+1), h5,-1(l=-1), h5,+2(l=+2), h5,-2(l=-2), h5,+3(l=+3), h5,-3(l=-3), h5,+4(l=+4), h5,-4(l=-4), h5,+5(l=+5), h5,-5(l=-5)" - ) - logger.debug( - "(makefun notation): h5,0(l=0), h5,+1(l=+1), h5,-1(l=-1), h5,+2(l=+2), h5,-2(l=-2), h5,+3(l=+3), h5,-3(l=-3), h5,+4(l=+4), h5,-4(l=-4), h5,+5(l=+5), h5,-5(l=-5)" - ) - reorder_index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] - reorder_m_list = [ - 0, - +1, - -1, - +2, - -2, - +3, - -3, - +4, - -4, - +5, - -5, - ] - reorder_l_list = [5] * 11 + elif current_ang_mom == 5: # h shell - elif current_ang_mom == 6: # i shell + if ao_cartesian != 0: + logger.debug( + "Cartesian notation for h shell is not implemented yet! Sorry." + ) + raise NotImplementedError - if ao_cartesian != 0: - logger.debug( - "Cartesian notation for i shell is not implemented yet! Sorry." - ) - raise NotImplementedError + else: + logger.debug("h shell/no permutation is needed.") + logger.debug( + "(trexio notation): h5,0(l=0), h5,+1(l=+1), h5,-1(l=-1), h5,+2(l=+2), h5,-2(l=-2), h5,+3(l=+3), h5,-3(l=-3), h5,+4(l=+4), h5,-4(l=-4), h5,+5(l=+5), h5,-5(l=-5)" + ) + logger.debug( + "(makefun notation): h5,0(l=0), h5,+1(l=+1), h5,-1(l=-1), h5,+2(l=+2), h5,-2(l=-2), h5,+3(l=+3), h5,-3(l=-3), h5,+4(l=+4), h5,-4(l=-4), h5,+5(l=+5), h5,-5(l=-5)" + ) + reorder_index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] + reorder_m_list = [ + 0, + +1, + -1, + +2, + -2, + +3, + -3, + +4, + -4, + +5, + -5, + ] + reorder_l_list = [5] * 11 + + elif current_ang_mom == 6: # i shell + + if ao_cartesian != 0: + logger.debug( + "Cartesian notation for i shell is not implemented yet! Sorry." + ) + raise NotImplementedError + + else: + logger.debug("i shell/no permutation is needed.") + logger.debug( + "(trexio notation): i6,0(l=0), i6,+1(l=+1), i6,-1(l=-1), i6,+2(l=+2), i6,-2(l=-2), i6,+3(l=+3), i6,-3(l=-3), i6,+4(l=+4), i6,-4(l=-4), i6,+5(l=+5), i6,-5(l=-5), i6,+6(l=+6), i6,-6(l=-6)" + ) + logger.debug( + "(makefun notation): i6,0(l=0), i6,+1(l=+1), i6,-1(l=-1), i6,+2(l=+2), i6,-2(l=-2), i6,+3(l=+3), i6,-3(l=-3), i6,+4(l=+4), i6,-4(l=-4), i6,+5(l=+5), i6,-5(l=-5), i6,+6(l=+6), i6,-6(l=-6)" + ) + reorder_index = [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + ] + reorder_m_list = [ + 0, + +1, + -1, + +2, + -2, + +3, + -3, + +4, + -4, + +5, + -5, + +6, + -6, + ] + reorder_l_list = [6] * 13 else: - logger.debug("i shell/no permutation is needed.") - logger.debug( - "(trexio notation): i6,0(l=0), i6,+1(l=+1), i6,-1(l=-1), i6,+2(l=+2), i6,-2(l=-2), i6,+3(l=+3), i6,-3(l=-3), i6,+4(l=+4), i6,-4(l=-4), i6,+5(l=+5), i6,-5(l=-5), i6,+6(l=+6), i6,-6(l=-6)" - ) - logger.debug( - "(makefun notation): i6,0(l=0), i6,+1(l=+1), i6,-1(l=-1), i6,+2(l=+2), i6,-2(l=-2), i6,+3(l=+3), i6,-3(l=-3), i6,+4(l=+4), i6,-4(l=-4), i6,+5(l=+5), i6,-5(l=-5), i6,+6(l=+6), i6,-6(l=-6)" + logger.error( + f" Angular momentum={ang_mom} is not implemented in TurboRVB!!!" ) - reorder_index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] - reorder_m_list = [ - 0, - +1, - -1, - +2, - -2, - +3, - -3, - +4, - -4, - +5, - -5, - +6, - -6, - ] - reorder_l_list = [6] * 13 + raise NotImplementedError - else: - logger.error( - f" Angular momentum={ang_mom} is not implemented in TurboRVB!!!" - ) - raise NotImplementedError + basis_index_list = [ + n + for n, v in enumerate(basis_shell_index) + if v == shell_index + ] - basis_index_list = [ - n for n, v in enumerate(basis_shell_index) if v == shell_index - ] + mo_coeff_local_buffer = [] + mo_exponent_local_buffer = [] + mo_nucleus_index_local_buffer = [] - mo_coeff_local_buffer = [] - mo_exponent_local_buffer = [] - mo_nucleus_index_local_buffer = [] + logger.debug(basis_index_list) + logger.debug(basis_nucleus_index) + for basis_index in basis_index_list: + mo_coeff_local_buffer.append( + mo_coefficient[mo_i_shifted][ao_i] + * basis_coefficient[basis_index] + ) + mo_exponent_local_buffer.append( + basis_exponent[basis_index] + ) + mo_nucleus_index_local_buffer.append( + basis_nucleus_index[shell_index] + ) - logger.debug(basis_index_list) - logger.debug(basis_nucleus_index) - for basis_index in basis_index_list: - mo_coeff_local_buffer.append( - mo_coefficient[mo_i][ao_i] * basis_coefficient[basis_index] + mo_coefficient_list_for_reordering.append( + mo_coeff_local_buffer ) - mo_exponent_local_buffer.append(basis_exponent[basis_index]) - mo_nucleus_index_local_buffer.append( - basis_nucleus_index[shell_index] + mo_exponent_list_for_reordering.append( + mo_exponent_local_buffer ) - - mo_coefficient_list_for_reordering.append(mo_coeff_local_buffer) - mo_exponent_list_for_reordering.append(mo_exponent_local_buffer) - mo_nucleus_index_list_for_reordering.append( - mo_nucleus_index_local_buffer - ) - - # write MOs!! - if len(mo_coefficient_list_for_reordering) == multiplicity: - logger.debug("--write MOs!!--") - mo_coeff_list_reordered = [ - mo_coefficient_list_for_reordering[i] - for i in reorder_index - ] - mo_exp_list_reordered = [ - mo_exponent_list_for_reordering[i] for i in reorder_index - ] - mo_nucleus_index_list_reordered = [ - mo_nucleus_index_list_for_reordering[i] - for i in reorder_index - ] - mo_m_list_reordered = reorder_m_list * len(basis_index_list) - mo_l_list_reordered = reorder_l_list * len(basis_index_list) - - # reordered again is needed here for contracted shell! - # because order is not px ... py ... pz but px,py,pz,px,py,pz ... - for ii in range(len(mo_coeff_list_reordered[0])): - for mo_coeff_reordered in mo_coeff_list_reordered: - mo_coefficient_list.append(mo_coeff_reordered[ii]) - for ii in range(len(mo_exp_list_reordered[0])): - for mo_exp_reordered in mo_exp_list_reordered: - mo_exponent_list.append(mo_exp_reordered[ii]) - for ii in range(len(mo_nucleus_index_list_reordered[0])): - for ( - mo_nucleus_index_reordered - ) in mo_nucleus_index_list_reordered: - mo_nucleus_index_list.append( - mo_nucleus_index_reordered[ii] - ) - """ - for ii in range(len(mo_m_list_reordered[0])): - for mo_m_reordered in mo_m_list_reordered: - mo_m_list.append(mo_m_reordered[ii]) - for ii in range(len(mo_angmom_list_reordered[0])): - for mo_angmom_reordered in mo_angmom_list_reordered: - mo_l_list.append(mo_angmom_reordered[ii]) - """ - mo_m_list += mo_m_list_reordered - mo_l_list += mo_l_list_reordered - - # reset buffer after writing MOs - mo_coefficient_list_for_reordering = [] - mo_exponent_list_for_reordering = [] - mo_nucleus_index_list_for_reordering = [] - - logger.debug(len(mo_exponent_list)) - logger.debug(len(mo_coefficient_list)) - logger.debug(len(mo_nucleus_index_list)) - logger.debug(len(mo_m_list)) - logger.debug(len(mo_l_list)) - - assert len(mo_coefficient_list) == len(mo_exponent_list) - assert len(mo_coefficient_list) == len(mo_nucleus_index_list) - assert len(mo_coefficient_list) == len(mo_m_list) - assert len(mo_coefficient_list) == len(mo_l_list) - - logger.debug("========================================") - logger.debug("End buffering MOs") - logger.debug("========================================") - - # remove duplicated exponents! - # Note: TurboRVB internally removes duplicated exponents. - logger.debug(mo_nucleus_index_list) - logger.debug(len(mo_nucleus_index_list)) - - ref_exponent_index_list = [] - removed_exponent_index_list = [] - for unique_mo_nucleus_index in set(mo_nucleus_index_list): - logger.debug( - f"unique_mo_nucleus_index = {unique_mo_nucleus_index}." - ) - mo_nucleus_index = [ - i - for i, x in enumerate(mo_nucleus_index_list) - if x == unique_mo_nucleus_index - ] - logger.debug(mo_nucleus_index) - mo_exponent_list_nucleus = [ - mo_exponent_list[i] for i in mo_nucleus_index - ] - mo_m_list_nucleus = [mo_m_list[i] for i in mo_nucleus_index] - mo_l_list_nucleus = [mo_l_list[i] for i in mo_nucleus_index] - logger.debug(mo_exponent_list_nucleus) - - unique_exp_list = list(set(mo_exponent_list_nucleus)) - logger.debug(mo_exponent_list_nucleus) - logger.debug(unique_exp_list) - - logger.debug(len(mo_exponent_list_nucleus)) - logger.debug(len(unique_exp_list)) - if len(unique_exp_list) != len(mo_exponent_list_nucleus): - logger.debug( - f"Duplicated exponents exist for nucleus index = {unique_mo_nucleus_index}." + mo_nucleus_index_list_for_reordering.append( + mo_nucleus_index_local_buffer ) - for unique_exp in unique_exp_list: - unique_exponent_index = [ - i - for i, x in enumerate(mo_exponent_list_nucleus) - if x == unique_exp - ] - mo_e_list_u = [ - mo_exponent_list_nucleus[i] - for i in unique_exponent_index + + # store MOs!! + if len(mo_coefficient_list_for_reordering) == multiplicity: + logger.debug("--write MOs!!--") + mo_coeff_list_reordered = [ + mo_coefficient_list_for_reordering[i] + for i in reorder_index ] - assert len(set(mo_e_list_u)) == 1 - mo_l_list_u = [ - mo_l_list_nucleus[i] for i in unique_exponent_index + mo_exp_list_reordered = [ + mo_exponent_list_for_reordering[i] + for i in reorder_index ] - mo_m_list_u = [ - mo_m_list_nucleus[i] for i in unique_exponent_index + mo_nucleus_index_list_reordered = [ + mo_nucleus_index_list_for_reordering[i] + for i in reorder_index ] + mo_m_list_reordered = reorder_m_list * len( + basis_index_list + ) + mo_l_list_reordered = reorder_l_list * len( + basis_index_list + ) - while ( - len(unique_exponent_index) > 0 - and len(mo_l_list_u) > 0 - and len(mo_m_list_u) > 0 - ): - ref_flag = False - removed_uu_list = [] - removed_exponent_index_list_ = [] - logger.debug(unique_exponent_index) - logger.debug(mo_l_list_nucleus) - logger.debug(mo_m_list_nucleus) - for uu, (e_i, mo_l, mo_m) in enumerate( - zip( - unique_exponent_index, mo_l_list_u, mo_m_list_u + # reordered again is needed here for contracted shell! + # because order is not px ... py ... pz but px,py,pz,px,py,pz ... + for ii in range(len(mo_coeff_list_reordered[0])): + for mo_coeff_reordered in mo_coeff_list_reordered: + mo_coefficient_list.append(mo_coeff_reordered[ii]) + for ii in range(len(mo_exp_list_reordered[0])): + for mo_exp_reordered in mo_exp_list_reordered: + mo_exponent_list.append(mo_exp_reordered[ii]) + for ii in range(len(mo_nucleus_index_list_reordered[0])): + for ( + mo_nucleus_index_reordered + ) in mo_nucleus_index_list_reordered: + mo_nucleus_index_list.append( + mo_nucleus_index_reordered[ii] ) - ): - if not ref_flag: - ref_e_i = e_i - mo_l_ref = mo_l - mo_m_ref = mo_m - removed_uu_list.append(uu) - ref_flag = True - else: - if mo_l_ref == mo_l and mo_m_ref == mo_m: - removed_uu_list.append(uu) - removed_exponent_index_list_.append(e_i) - - logger.debug(removed_uu_list) - for kk in reversed(sorted(removed_uu_list)): - unique_exponent_index.pop(kk) - mo_l_list_u.pop(kk) - mo_m_list_u.pop(kk) + """ + for ii in range(len(mo_m_list_reordered[0])): + for mo_m_reordered in mo_m_list_reordered: + mo_m_list.append(mo_m_reordered[ii]) + for ii in range(len(mo_angmom_list_reordered[0])): + for mo_angmom_reordered in mo_angmom_list_reordered: + mo_l_list.append(mo_angmom_reordered[ii]) + """ + mo_m_list += mo_m_list_reordered + mo_l_list += mo_l_list_reordered + + # reset buffer after writing MOs + mo_coefficient_list_for_reordering = [] + mo_exponent_list_for_reordering = [] + mo_nucleus_index_list_for_reordering = [] + + assert len(mo_coefficient_list) == len(mo_exponent_list) + assert len(mo_coefficient_list) == len(mo_nucleus_index_list) + assert len(mo_coefficient_list) == len(mo_m_list) + assert len(mo_coefficient_list) == len(mo_l_list) + + # ======================================== + # End buffering MOs + # ======================================== + + # remove duplicated exponents! + # Note: TurboRVB internally removes duplicated exponents. + + ref_exponent_index_list = [] + removed_exponent_index_list = [] + for unique_mo_nucleus_index in set(mo_nucleus_index_list): + mo_nucleus_index = [ + i + for i, x in enumerate(mo_nucleus_index_list) + if x == unique_mo_nucleus_index + ] + mo_exponent_list_nucleus = [ + mo_exponent_list[i] for i in mo_nucleus_index + ] + mo_m_list_nucleus = [mo_m_list[i] for i in mo_nucleus_index] + mo_l_list_nucleus = [mo_l_list[i] for i in mo_nucleus_index] - logger.debug(ref_e_i) - logger.debug(removed_exponent_index_list_) - logger.debug(mo_nucleus_index[ref_e_i]) - logger.debug( - [ - mo_nucleus_index[pp] - for pp in removed_exponent_index_list_ - ] - ) - ref_exponent_index_list.append( - mo_nucleus_index[ref_e_i] - ) - removed_exponent_index_list.append( - [ - mo_nucleus_index[pp] - for pp in removed_exponent_index_list_ - ] - ) + unique_exp_list = list(set(mo_exponent_list_nucleus)) - else: - logger.debug( - f"No duplicated exponent is found for nucleus index ={unique_mo_nucleus_index}" - ) + if len(unique_exp_list) != len(mo_exponent_list_nucleus): + for unique_exp in unique_exp_list: + unique_exponent_index = [ + i + for i, x in enumerate(mo_exponent_list_nucleus) + if x == unique_exp + ] + mo_e_list_u = [ + mo_exponent_list_nucleus[i] + for i in unique_exponent_index + ] + assert len(set(mo_e_list_u)) == 1 + mo_l_list_u = [ + mo_l_list_nucleus[i] for i in unique_exponent_index + ] + mo_m_list_u = [ + mo_m_list_nucleus[i] for i in unique_exponent_index + ] - # Finally, here, duplicated exponents are removed. - # because of using pop here. - # if pop is used inside the above for loop, indices change during the for loop. - pop_index_list = [] - for ref_exponent_index, removed_exponent_index in zip( - ref_exponent_index_list, removed_exponent_index_list - ): - for r_index in reversed(sorted(removed_exponent_index)): - mo_coefficient_list[ref_exponent_index] += mo_coefficient_list[ - r_index - ] - pop_index_list.append(r_index) - logger.debug(sorted(pop_index_list)) - logger.debug(len(sorted(pop_index_list))) - for pop_index in reversed(sorted(pop_index_list)): - mo_coefficient_list.pop(pop_index) + while ( + len(unique_exponent_index) > 0 + and len(mo_l_list_u) > 0 + and len(mo_m_list_u) > 0 + ): + ref_flag = False + removed_uu_list = [] + removed_exponent_index_list_ = [] + logger.debug(unique_exponent_index) + logger.debug(mo_l_list_nucleus) + logger.debug(mo_m_list_nucleus) + for uu, (e_i, mo_l, mo_m) in enumerate( + zip( + unique_exponent_index, + mo_l_list_u, + mo_m_list_u, + ) + ): + if not ref_flag: + ref_e_i = e_i + mo_l_ref = mo_l + mo_m_ref = mo_m + removed_uu_list.append(uu) + ref_flag = True + else: + if mo_l_ref == mo_l and mo_m_ref == mo_m: + removed_uu_list.append(uu) + removed_exponent_index_list_.append( + e_i + ) + + logger.debug(removed_uu_list) + for kk in reversed(sorted(removed_uu_list)): + unique_exponent_index.pop(kk) + mo_l_list_u.pop(kk) + mo_m_list_u.pop(kk) + + ref_exponent_index_list.append( + mo_nucleus_index[ref_e_i] + ) + removed_exponent_index_list.append( + [ + mo_nucleus_index[pp] + for pp in removed_exponent_index_list_ + ] + ) - mo_coefficient_turbo.append(mo_coefficient_list) - mo_exponent_turbo.append(mo_exponent_list) + else: + logger.debug( + f"No duplicated exponent is found for nucleus index ={unique_mo_nucleus_index}" + ) - # fort.10 MO replace - logger.debug(len(io_fort10.f10detbasissets.mo_coefficient[0])) - logger.debug(len(mo_coefficient_turbo[0])) + # Finally, here, duplicated exponents are removed. + # because of using pop here. + # if pop is used inside the above for loop, indices change during the for loop. + pop_index_list = [] + for ref_exponent_index, removed_exponent_index in zip( + ref_exponent_index_list, removed_exponent_index_list + ): + for r_index in reversed(sorted(removed_exponent_index)): + mo_coefficient_list[ + ref_exponent_index + ] += mo_coefficient_list[r_index] + pop_index_list.append(r_index) + + for pop_index in reversed(sorted(pop_index_list)): + mo_coefficient_list.pop(pop_index) + + mo_coefficient_turbo.append(mo_coefficient_list) + mo_exponent_turbo.append(mo_exponent_list) + + # mo_exponent_turbo is duplicated in the spin-unrestricted case + # remove the duplication here + if not spin_restricted: + if ( + mo_exponent_turbo[0:mo_num_use] + != mo_exponent_turbo[mo_num_use : 2 * mo_num_use] + ): + logger.error( + "mo_exponent_turbo is not consistent between alpha and beta spins." + ) + raise ValueError + mo_exponent_turbo = mo_exponent_turbo[0:mo_num_use] # molecular orbital swapped, spin polarized cases. - mo_coefficient_turbo_unpaired = [] - pop_lst = [ - num_ele_dn + nel_diff for nel_diff in range(0, num_ele_up - num_ele_dn) - ] - for p in reversed(pop_lst): - logger.debug("molecular orbital swapped") - mo_coefficient_turbo_unpaired.append(mo_coefficient_turbo.pop(p)) - for m in reversed(mo_coefficient_turbo_unpaired): - mo_coefficient_turbo.append(m) + if spin_restricted: + logger.info("Molecular orbitals are swapped (spin-resticted case).") + # spin-restricted case, here, the MOs are reordered such that + # [a,a,a,a(unpaired),a(unpaired)....a] + # -> + # # [a,a,a... (paired),a,a (unpaired)] + mo_coefficient_turbo_unpaired = [] + pop_lst = [ + num_ele_dn + nel_diff + for nel_diff in range(0, num_ele_up - num_ele_dn) + ] + for p in reversed(pop_lst): + mo_coefficient_turbo_unpaired.append(mo_coefficient_turbo.pop(p)) + for m in reversed(mo_coefficient_turbo_unpaired): + mo_coefficient_turbo.append(m) + else: + logger.info("Molecular orbitals are swapped (spin-unresticted case).") + # spin-unrestricted case, here, the MOs are reordered such that + # [a,a,a,a,a....a, b,b,b,b,b,b,b,....b] + # -> + # [b,a,b,a,b,a,b......a (paired),a,a,a,a (unpaired)] + # molecular orbital swapped, spin polarized cases. + mo_coefficient_turbo_unpaired = [] + # alpha spins + pop_lst = [ + num_ele_dn + nel_diff + for nel_diff in range(0, num_ele_up - num_ele_dn) + ] + + for p in reversed(pop_lst): + mo_coefficient_turbo_unpaired.append(mo_coefficient_turbo.pop(p)) + # beta spins + # note!! + # there are two choices which MOs are removed in the beta spins + # 1. the last 'nel_diff' beta orbitals are removed. + # 2. the beta MOs corresponding to unpaired alpha MOs are removed. + #for _ in range(0, num_ele_up - num_ele_dn): # this is the choice 1. + # mo_coefficient_turbo.pop(-1) + for _ in range(0, num_ele_up - num_ele_dn): # this is the choice 2. + mo_coefficient_turbo.pop(mo_num_use - len(pop_lst) + num_ele_dn) + # reordering alpha and beta spins + mo_coefficient_turbo_new = [] + for i in range(0, int(len(mo_coefficient_turbo) / 2)): + mo_coefficient_turbo_new.append( + mo_coefficient_turbo[i + int(len(mo_coefficient_turbo) / 2)] + ) # beta spin + mo_coefficient_turbo_new.append( + mo_coefficient_turbo[i + 0] + ) # alpha spin + mo_coefficient_turbo = mo_coefficient_turbo_new + for m in reversed(mo_coefficient_turbo_unpaired): + mo_coefficient_turbo.append(m) # fort.10 MO replace logger.info("Writing obtained MOs to fort.10....(It might take a while).") if not complex_flag: + # logger.info(len(mo_coefficient_turbo)) + # logger.info(len(io_fort10.f10detbasissets.mo_coefficient)) io_fort10.f10detbasissets.mo_coefficient = mo_coefficient_turbo else: # complex case, # mo_coefficient_turbo -> mo_coefficient_turbo_real, mo_coefficient_turbo_imag mo_coefficient_turbo_real = [] mo_coefficient_turbo_imag = [] - for mo__ in mo_coefficient_turbo: + for mo_counter, mo__ in enumerate(mo_coefficient_turbo): mo_real_b = [] mo_imag_b = [] for coeff in mo__: mo_real_b.append(coeff.real) mo_imag_b.append(coeff.imag) - - mo_real_b_up = list(np.array(mo_real_b) * +1) - mo_imag_b_up = list(np.array(mo_imag_b) * +1) # up phase is + - # these commented lines are wrong!! In general, the wf does not symmetric with respect to the time reversal except for TRIM points. - # mo_real_b_dn=list(np.array(mo_real_b) * +1) # real part is the same as the up spin. - # mo_imag_b_dn=list(np.array(mo_imag_b) * -1) # dn phase is -. because the opposite phase is attached in turbo with option double k-grid=.true. - mo_real_b_dn = list(np.array(mo_real_b) * +1) - mo_imag_b_dn = list(np.array(mo_imag_b) * +1) - - # turbo seems to treat up and dn MOs individually for complex cases!! - # trexio currently supports only restricted open or closed shells, => This is a TODO - # so, here just duplicate the MOs. => is it ok?? - # Yes. the order seems correct, i.e, up, dn, up, dn.... - # but, let me see,,, for spin-polarized cases?? - if num_ele_up != num_ele_dn: - logger.error( - "spin-polarized case for a complex system is not implemented yet." - ) - raise NotImplementedError - - mo_coefficient_turbo_real.append(mo_real_b_up) # up - mo_coefficient_turbo_imag.append(mo_imag_b_up) # up - mo_coefficient_turbo_real.append(mo_real_b_dn) # dn - mo_coefficient_turbo_imag.append(mo_imag_b_dn) # dn + + if spin_restricted: + mo_real_b_up = list(np.array(mo_real_b) * +1) + mo_imag_b_up = list(np.array(mo_imag_b) * +1) # up phase is + + # these commented lines are wrong!! In general, the wf does not symmetric with respect to the time reversal except for TRIM points. + # mo_real_b_dn=list(np.array(mo_real_b) * +1) # real part is the same as the up spin. + # mo_imag_b_dn=list(np.array(mo_imag_b) * -1) # dn phase is -. because the opposite phase is attached in turbo with option double k-grid=.true. + mo_real_b_dn = list(np.array(mo_real_b) * +1) + mo_imag_b_dn = list(np.array(mo_imag_b) * +1) + + mo_coefficient_turbo_real.append(mo_real_b_dn) # dn + mo_coefficient_turbo_imag.append(mo_imag_b_dn) # dn + # because mo_dn is not needed for unpaired MOs. + # if num_ele_up - num_ele_dn == 0: the following condition is always true. + if (len(mo_coefficient_turbo) - (num_ele_up - num_ele_dn)) >= mo_counter + 1: + mo_coefficient_turbo_real.append(mo_real_b_up) # up + mo_coefficient_turbo_imag.append(mo_imag_b_up) # up + else: # if spin_restricted==False + mo_coefficient_turbo_real.append(mo_real_b) + mo_coefficient_turbo_imag.append(mo_imag_b) logger.debug(mo_coefficient_turbo) - logger.debug(len(io_fort10.f10detbasissets.mo_coefficient)) - logger.debug(len(mo_coefficient_turbo_real)) - logger.debug(len(io_fort10.f10detbasissets.mo_coefficient_imag)) - logger.debug(len(mo_coefficient_turbo_imag)) + logger.info(f"fort10mo_real={len(io_fort10.f10detbasissets.mo_coefficient)}") + logger.info(f"trexmo_real={len(mo_coefficient_turbo_real)}") + logger.info(f"fort10mo_real={len(io_fort10.f10detbasissets.mo_coefficient_imag)}") + logger.info(f"trexmo_real={len(mo_coefficient_turbo_imag)}") io_fort10.f10detbasissets.mo_coefficient = mo_coefficient_turbo_real io_fort10.f10detbasissets.mo_coefficient_imag = ( mo_coefficient_turbo_imag @@ -1101,7 +1330,7 @@ def checker(choice): trexio_to_turborvb_wf( trexio_file=args.trexio_file, cleanup=args.cleanup, - max_occ_conv=0.0, + max_occ_conv=0.01, jas_basis_sets=jas_basis_sets, ) diff --git a/turbogenius/trexio_wrapper.py b/turbogenius/trexio_wrapper.py index 996c915..0f8f6bd 100644 --- a/turbogenius/trexio_wrapper.py +++ b/turbogenius/trexio_wrapper.py @@ -32,7 +32,7 @@ class Trexio_wrapper_r: """ - def __init__(self, trexio_file): + def __init__(self, trexio_file: str): # prefix and file names logger.info(f"TREXIO file = {trexio_file}") @@ -83,16 +83,17 @@ def __init__(self, trexio_file): self.basis_prim_factor = trexio.read_basis_prim_factor(file_r) # Pseudo potentials info - self.ecp_max_ang_mom_plus_1 = trexio.read_ecp_max_ang_mom_plus_1( - file_r - ) - self.ecp_z_core = trexio.read_ecp_z_core(file_r) - self.ecp_num = trexio.read_ecp_num(file_r) - self.ecp_ang_mom = trexio.read_ecp_ang_mom(file_r) - self.ecp_nucleus_index = trexio.read_ecp_nucleus_index(file_r) - self.ecp_exponent = trexio.read_ecp_exponent(file_r) - self.ecp_coefficient = trexio.read_ecp_coefficient(file_r) - self.ecp_power = trexio.read_ecp_power(file_r) + if trexio.has_ecp_num(file_r): + self.ecp_max_ang_mom_plus_1 = trexio.read_ecp_max_ang_mom_plus_1( + file_r + ) + self.ecp_z_core = trexio.read_ecp_z_core(file_r) + self.ecp_num = trexio.read_ecp_num(file_r) + self.ecp_ang_mom = trexio.read_ecp_ang_mom(file_r) + self.ecp_nucleus_index = trexio.read_ecp_nucleus_index(file_r) + self.ecp_exponent = trexio.read_ecp_exponent(file_r) + self.ecp_coefficient = trexio.read_ecp_coefficient(file_r) + self.ecp_power = trexio.read_ecp_power(file_r) # ao info self.ao_cartesian = trexio.read_ao_cartesian(file_r) @@ -105,6 +106,10 @@ def __init__(self, trexio_file): self.mo_num = trexio.read_mo_num(file_r) self.mo_occupation = trexio.read_mo_occupation(file_r) self.mo_coefficient = trexio.read_mo_coefficient(file_r) + try: + self.mo_spin = trexio.read_mo_spin(file_r) + except: # backward compatibility + self.mo_spin = [0 for _ in range(self.mo_num)] if trexio.has_mo_coefficient_im(file_r): logger.info("The WF is complex") self.mo_coefficient_imag = trexio.read_mo_coefficient_im(file_r) @@ -115,35 +120,6 @@ def __init__(self, trexio_file): file_r.close() - """ - def write_to_turbowf( - self, - jas_basis_sets=Jas_Basis_sets(), - max_occ_conv: int = 0, - mo_num_conv: int = -1, - only_mol: bool = True, - cleanup: bool = True, - ) -> None: - - Convert trexio file to TurboRVB WF file (fort.10) - - Args: - jas_basis_sets (Jas_basis_sets): Jastrow basis sets added to the TREXIO WF. - max_occ_conv (int): maximum occ used for the conv, not used with mo_num - mo_num_conv (int): num mo used for the conv, not used with max occ - only_mol (bool): if True, only moleculer orbitals option = True in convertfort10mol - cleanup (bool): clean up temporary files - - trexio_to_turborvb_wf( - trexio_file=self.trexio_file, - jas_basis_sets=jas_basis_sets, - max_occ_conv=max_occ_conv, - mo_num_conv=max_occ_conv, - only_mol=only_mol, - cleanup=cleanup, - ) - """ - if __name__ == "__main__": import sys diff --git a/turbogenius/turbo_genius_cli.py b/turbogenius/turbo_genius_cli.py index 14ca566..0c74731 100755 --- a/turbogenius/turbo_genius_cli.py +++ b/turbogenius/turbo_genius_cli.py @@ -891,8 +891,7 @@ def vmcopt( else: logger.info("Job was failure. See the output file.") return - if plot_graph: - vmcopt_genius.plot_energy_and_devmax(interactive=plot_interactive) + vmcopt_genius.plot_energy_and_devmax(interactive=plot_interactive) if plot_graph: vmcopt_genius.plot_parameters_history(interactive=plot_interactive) vmcopt_genius.average( @@ -1577,7 +1576,11 @@ def convertpfaff( @cli.command(short_help="convert wavefunction") @decorate_grpost @click.option( - "-to", "to_ansatz", help="Specify to ansatz", default="agps", type=str + "-to", + "to_ansatz", + help="Specify to ansatz", + default="agps", + type=click.Choice(["sd", "agps", "agpu", "pf"], case_sensitive=False), ) @click.option( "-hyb", @@ -1622,7 +1625,8 @@ def convertwf( os.chdir(root_dir) number_of_additional_hybrid_orbitals = list(map(int, hybrid_orbitals)) - wavefunction = Wavefunction(fort10="fort.10") + wavefunction = Wavefunction() + wavefunction.read_from_fort10(fort10="fort.10") if to_ansatz == "agps": logger.info("convert the wf to agps") @@ -1630,7 +1634,7 @@ def convertwf( grid_size=grid_size, additional_hyb=number_of_additional_hybrid_orbitals, nosym=nosymmetry, - clean_flag=False, + clean_flag=True, ) elif to_ansatz == "agpu": @@ -1639,7 +1643,7 @@ def convertwf( grid_size=grid_size, additional_hyb=number_of_additional_hybrid_orbitals, nosym=nosymmetry, - clean_flag=False, + clean_flag=True, ) elif to_ansatz == "sd": @@ -1653,7 +1657,7 @@ def convertwf( grid_size=grid_size, rotate_angle=rotate_angle, nosym=nosymmetry, - clean_flag=False, + clean_flag=True, ) diff --git a/turbogenius/utils_workflows/utility.py b/turbogenius/utils_workflows/utility.py index c538447..8b3850a 100755 --- a/turbogenius/utils_workflows/utility.py +++ b/turbogenius/utils_workflows/utility.py @@ -17,7 +17,7 @@ def prompt(text, checker): return output -def get_nonlocalmoves_setting(nonlocalmoves): +def get_nonlocalmoves_setting(nonlocalmoves: str): if nonlocalmoves == "tmove": typereg = 0 npow = 0.0 @@ -38,16 +38,16 @@ def get_nonlocalmoves_setting(nonlocalmoves): def get_optimizer_flags( - optimizer="sr", - opt_onebody=True, - opt_twobody=True, - opt_det_mat=True, - opt_jas_mat=True, - opt_det_basis_exp=True, - opt_jas_basis_exp=True, - opt_det_basis_coeff=True, - opt_jas_basis_coeff=True, - qmc_type="vmc", + optimizer: str = "sr", + opt_onebody: bool = True, + opt_twobody: bool = True, + opt_det_mat: bool = True, + opt_jas_mat: bool = True, + opt_det_basis_exp: bool = True, + opt_jas_basis_exp: bool = True, + opt_det_basis_coeff: bool = True, + opt_jas_basis_coeff: bool = True, + qmc_type: str = "vmc", ): logger.debug(optimizer) @@ -176,7 +176,7 @@ def get_optimizer_flags( def return_optimizer_number( - optimizer="sr", qmc_type="vmc", opt_exponent=False + optimizer: str = "sr", qmc_type: str = "vmc", opt_exponent: bool = False ): if optimizer == "lr": if qmc_type == "vmc": diff --git a/turbogenius/vmc_genius.py b/turbogenius/vmc_genius.py index 8c757dc..0cdb315 100755 --- a/turbogenius/vmc_genius.py +++ b/turbogenius/vmc_genius.py @@ -12,6 +12,7 @@ # python modules import os +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -46,10 +47,13 @@ def __init__( num_walkers: int = -1, # default -1 -> num of MPI process. maxtime: int = 172800, twist_average: bool = False, - kpoints: list = [1, 1, 1, 0, 0, 0], + kpoints: Optional[list] = None, force_calc_flag: bool = True, ): + if kpoints is None: + kpoints = [1, 1, 1, 0, 0, 0] + self.force_calc_flag = force_calc_flag self.vmcsteps = vmcsteps self.num_walkers = num_walkers @@ -219,7 +223,7 @@ def store_result( self, bin_block: int = 10, warmupblocks: int = 5, - output_names: list = ["out_vmc"], + output_names: Optional[list] = None, rerun: bool = False, ) -> bool: """ @@ -232,6 +236,8 @@ def store_result( output_names (list): a list of output file names rerun (bool): if true, compute energy and force again even if there are energy and force files. """ + if output_names is None: + output_names = ["out_vmc"] self.estimated_time_for_1_generation = ( self.get_estimated_time_for_1_generation(output_names=output_names) ) @@ -259,7 +265,7 @@ def compute_energy_and_forces( ) def get_estimated_time_for_1_generation( - self, output_names: list = ["out_vmc"] + self, output_names: Optional[list] = None ) -> float: """ This procedure stores estimated_time_for_1_generation. @@ -270,11 +276,13 @@ def get_estimated_time_for_1_generation( Return: float: estimated_time_for_1_generation. """ + if output_names is None: + output_names = ["out_vmc"] return self.vmc.get_estimated_time_for_1_generation( output_names=output_names ) - def check_results(self, output_names: list = ["out_vmc"]) -> bool: + def check_results(self, output_names: Optional[list] = None) -> bool: """ Check the result. @@ -283,6 +291,8 @@ def check_results(self, output_names: list = ["out_vmc"]) -> bool: Return: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_vmc"] return self.vmc.check_results(output_names=output_names) diff --git a/turbogenius/vmc_opt_genius.py b/turbogenius/vmc_opt_genius.py index 2bc26da..56f5aaa 100755 --- a/turbogenius/vmc_opt_genius.py +++ b/turbogenius/vmc_opt_genius.py @@ -12,6 +12,7 @@ # python modules import os +from typing import Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -80,9 +81,12 @@ def __init__( opt_structure: bool = False, str_learning_rate: float = 1.0e-6, twist_average: bool = False, - kpoints: list = [1, 1, 1, 0, 0, 0], + kpoints: Optional[list] = None, ): + if kpoints is None: + kpoints = [1, 1, 1, 0, 0, 0] + self.fort10 = fort10 self.twist_average = twist_average self.kpoints = kpoints @@ -364,7 +368,7 @@ def run( flags = self.vmcopt.check_results(output_names=[output_name]) assert all(flags) - def check_results(self, output_names: list = ["out_min"]) -> bool: + def check_results(self, output_names: Optional[list] = None) -> bool: """ Check the result. @@ -373,10 +377,12 @@ def check_results(self, output_names: list = ["out_min"]) -> bool: Returns: bool: True if all the runs were successful, False if an error is detected in the files. """ + if output_names is None: + output_names = ["out_min"] return self.vmcopt.check_results(output_names=output_names) def plot_energy_and_devmax( - self, output_names: list = ["out_min"], interactive: bool = True + self, output_names: Optional[list] = None, interactive: bool = True ) -> None: """ plot energy and devmax @@ -385,17 +391,21 @@ def plot_energy_and_devmax( output_names (list): a list of output file names interactive (bool): flag for an interactive plot """ + if output_names is None: + output_names = ["out_min"] self.vmcopt.plot_energy_and_devmax( output_names=output_names, interactive=interactive ) - def store_result(self, output_names: list = ["out_min"]) -> None: + def store_result(self, output_names: Optional[list] = None) -> None: """ Store results. energy, energy_error, and estimated_time_for_1_generation are stored in this class. Args: output_names (list): a list of output file names """ + if output_names is None: + output_names = ["out_min"] self.energy, self.energy_error = self.get_energy( output_names=output_names ) @@ -407,7 +417,7 @@ def average( self, optwarmupsteps: int, input_name: str = "datasmin.input", - output_names: list = ["out_min"], + output_names: Optional[list] = None, graph_plot: bool = False, ) -> None: """ @@ -419,6 +429,8 @@ def average( output_names (list): a list of output file names graph_plot (bool): Flag for plotting a graph """ + if output_names is None: + output_names = ["out_min"] current_dir = os.getcwd() os.chdir(current_dir) if self.twist_average: @@ -476,7 +488,7 @@ def average( os.chdir(current_dir) """ - def get_energy(self, output_names: list = ["out_min"]) -> list: + def get_energy(self, output_names: Optional[list] = None) -> list: """ return energy list @@ -486,10 +498,12 @@ def get_energy(self, output_names: list = ["out_min"]) -> list: Return: list: a list of history of energies. """ + if output_names is None: + output_names = ["out_min"] return self.vmcopt.get_energy(output_names=output_names) def get_estimated_time_for_1_generation( - self, output_names: list = ["out_min"] + self, output_names: Optional[list] = None ) -> float: """ This procedure stores estimated_time_for_1_generation. @@ -500,6 +514,8 @@ def get_estimated_time_for_1_generation( Return: float: estimated_time_for_1_generation. """ + if output_names is None: + output_names = ["out_min"] return self.vmcopt.get_estimated_time_for_1_generation( output_names=output_names ) diff --git a/turbogenius/wavefunction.py b/turbogenius/wavefunction.py index 694bafd..b565f26 100755 --- a/turbogenius/wavefunction.py +++ b/turbogenius/wavefunction.py @@ -16,7 +16,7 @@ import os import shutil import numpy as np -from typing import Union +from typing import Union, Optional # Logger from logging import getLogger, StreamHandler, Formatter @@ -309,7 +309,7 @@ def to_sd(self, grid_size: float = 0.10, clean_flag: bool = False) -> None: def to_agps( self, grid_size: float = 0.10, - additional_hyb: list = [], + additional_hyb: Optional[list] = None, nosym: float = False, clean_flag: float = False, ) -> None: @@ -323,13 +323,15 @@ def to_agps( clean_flag (bool): cleaning temporary files """ + if additional_hyb is None: + additional_hyb = [] if not self.read_flag: logger.warning( "WF file is not read yet! Please read a WF file first." ) return - self.io_fort10.to_agp( + self.to_agp( triplet=False, grid_size=grid_size, additional_hyb=additional_hyb, @@ -342,7 +344,7 @@ def to_agps( def to_agpu( self, grid_size: float = 0.10, - additional_hyb: list = [], + additional_hyb: Optional[list] = None, nosym: bool = False, clean_flag: bool = False, ) -> None: @@ -356,13 +358,15 @@ def to_agpu( clean_flag (bool): cleaning temporary files """ + if additional_hyb is None: + additional_hyb = [] if not self.read_flag: logger.warning( "WF file is not read yet! Please read a WF file first." ) return - self.io_fort10.to_agp( + self.to_agp( triplet=True, grid_size=grid_size, additional_hyb=additional_hyb, @@ -377,7 +381,7 @@ def to_agp( triplet: bool = False, pfaffian_flag: bool = False, grid_size: float = 0.10, - additional_hyb: list = [], + additional_hyb: Optional[list] = None, nosym: bool = False, clean_flag: bool = False, only_generate_template: bool = False, @@ -394,6 +398,8 @@ def to_agp( clean_flag (bool): cleaning temporary files """ + if additional_hyb is None: + additional_hyb = [] if not self.read_flag: logger.warning( "WF file is not read yet! Please read a WF file first." @@ -536,23 +542,25 @@ def to_agp( shutil.move("fort.10_new", "fort.10_out") if only_generate_template: - return - - # convertfort10 - convertfort10_genius = Convertfort10_genius( - in_fort10="fort.10_in", - out_fort10="fort.10_out", - grid_size=grid_size, - ) + logger.warning("A template AGP file, fort.10_out, is generated.") + else: + # convertfort10 + convertfort10_genius = Convertfort10_genius( + in_fort10="fort.10_in", + out_fort10="fort.10_out", + grid_size=grid_size, + ) - convertfort10_genius.generate_input(input_name="convertfort10.input") - convertfort10_genius.run( - input_name="convertfort10.input", output_name="out_conv" - ) - shutil.move(self.io_fort10.fort10, "fort.10_bak") - shutil.move("fort.10_new", "fort.10") - shutil.copy("fort.10_in", "fort.10_new") - copy_jastrow(fort10_to="fort.10", fort10_from="fort.10_new") + convertfort10_genius.generate_input( + input_name="convertfort10.input" + ) + convertfort10_genius.run( + input_name="convertfort10.input", output_name="out_conv" + ) + shutil.move(self.io_fort10.fort10, "fort.10_bak") + shutil.move("fort.10_new", "fort.10") + shutil.copy("fort.10_in", "fort.10_new") + copy_jastrow(fort10_to="fort.10", fort10_from="fort.10_new") if clean_flag: os.remove("fort.10_new")