diff --git a/.github/workflows/python-testsuite.yml b/.github/workflows/python-testsuite.yml index eebff323..cd0bbb82 100644 --- a/.github/workflows/python-testsuite.yml +++ b/.github/workflows/python-testsuite.yml @@ -13,11 +13,11 @@ jobs: build: runs-on: ubuntu-latest - + strategy: matrix: python-version: [3.8] - + services: rabbitmq: image: rabbitmq:latest @@ -38,6 +38,7 @@ jobs: pip install flake8 pytest~=6.0 pgtest~=1.3 aiida-core~=2.3 aiida-quantumespresso~=4.3 if [ ${{matrix.python-version}} -eq 2.7 ]; then pip install -r requirements2.txt; else pip install -r requirements.txt; fi + aiida-pseudo install - name: Lint with flake8 run: | @@ -54,7 +55,7 @@ jobs: cd CellConstructor python setup.py install --user cd .. - + python setup.py install --user # Install julia requirements diff --git a/Modules/AdvancedCalculations.py b/Modules/AdvancedCalculations.py index cf56ca8f..1b48c756 100644 --- a/Modules/AdvancedCalculations.py +++ b/Modules/AdvancedCalculations.py @@ -11,7 +11,7 @@ This module contains useful functions to do post-processing analysis of the SSCHA. For example, it implements custom cluster calculators to compute the electron-phonon effect on -absorbption and the bandgap. +absorption and the bandgap. ''' @@ -23,7 +23,7 @@ class OpticalQECluster(sscha.Cluster.Cluster): ''' def __init__(self, new_k_grid = None, random_offset = True, epsilon_data = None, - epsilon_binary = 'epsilon.x -npool NPOOL -i PREFIX.pwi > PREFIX.pwo', + epsilon_binary = 'epsilon.x -npool NPOOL -i PREFIX.pwi > PREFIX.pwo', **kwargs): ''' Initialize the cluster object. @@ -40,7 +40,7 @@ def __init__(self, new_k_grid = None, random_offset = True, epsilon_data = None, epsilon.x file epsilon_binary : string The path to the epsilon.x binary inside the cluster. - **kwargs : + **kwargs : All other arguments to be passed to the cluster. ''' @@ -67,7 +67,7 @@ def __init__(self, new_k_grid = None, random_offset = True, epsilon_data = None, 'wmax' : 40, 'nw' : 10000, 'temperature' : 300 - } + } } if epsilon_data is not None: self.epsilon_data = epsilon_data @@ -79,10 +79,10 @@ def __setattr__(self, __name, __value): super().__setattr__(__name, __value) # Always regenerate the kpts - if __name == 'new_k_grid' and not __value is None: + if __name == 'new_k_grid' and not __value is None: assert len(__value) == 3, 'Error, new_k_grid must be a tuple with 3 elements' self.generate_kpts() - + def generate_kpts(self): ''' @@ -119,7 +119,7 @@ def get_execution_command(self, label): # Get the MPI command replacing NPROC new_mpicmd = self.mpi_cmd.replace("NPROC", str(self.n_cpu)) - + # Replace the NPOOL variable and the PREFIX in the binary binary = self.binary.replace("NPOOL", str(self.n_pool)).replace("PREFIX", label) binary_nscf = self.binary.replace("NPOOL", str(self.n_pool)).replace("PREFIX", label + '_nscf') @@ -147,6 +147,18 @@ def get_execution_command(self, label): def read_results(self, calc, label): ''' Get the results + + Parameters + ---------- + calc : the ASE or CellConstructor calculator. + In this case, it works with quantum espresso + labels : List of strings + The unique name of this calculation + Returns + ------- + results : array + An array containing the real part of epsilon, + and the mean of the real and imaginary parst of epsilon. ''' results = super().read_results(calc, label) @@ -176,7 +188,31 @@ def read_results(self, calc, label): def prepare_input_file(self, structures, calc, labels): ''' Prepare the input files for the cluster - ''' + + This is specific for quantum espresso and must be inherit and replaced for + other calculators. + + This crates the input files in the local working directory + self.local_workdir and it returns the list of all the files generated. + + + Parameters + ---------- + structures : List of cellconstructor.Structure.Structure + The atomic structures. + calc : the ASE or CellConstructor calculator. + In this case, it works with quantum espresso + labels : List of strings + The unique name of this calculation + + Returns + ------- + List_of_input : list + List of strings containing all the input files + List_of_output : list + List of strings containing the output files expected + for the calculation + ''' # Prepare the input file @@ -185,7 +221,7 @@ def prepare_input_file(self, structures, calc, labels): for i, (label, structure) in enumerate(zip(labels, structures)): # Avoid thread conflict self.lock.acquire() - + try: calc.set_directory(self.local_workdir) PREFIX = label @@ -204,7 +240,7 @@ def prepare_input_file(self, structures, calc, labels): input_file = '{}.pwi'.format(label) output_file = '{}.pwo'.format(label) - list_of_inputs.append(input_file) + list_of_inputs.append(input_file) list_of_outputs.append(output_file) # prepare the nscf calculation @@ -218,7 +254,7 @@ def prepare_input_file(self, structures, calc, labels): calc.write_input(structure) input_file = '{}.pwi'.format(new_label) output_file = '{}.pwo'.format(new_label) - list_of_inputs.append(input_file) + list_of_inputs.append(input_file) list_of_outputs.append(output_file) @@ -239,8 +275,8 @@ def prepare_input_file(self, structures, calc, labels): print('fname: {} prepared'.format(eps_in_filename)) - - list_of_inputs.append(input_file) + + list_of_inputs.append(input_file) list_of_outputs.append(output_file) # Append also the imaginary and real part of epsilon and sigma @@ -263,11 +299,11 @@ def prepare_input_file(self, structures, calc, labels): print(e) # Release the lock on the threads - self.lock.release() + self.lock.release() print('THREAD: {} inputs: {} outputs: {}'.format(threading.get_native_id(), list_of_inputs, list_of_outputs)) - - + + return list_of_inputs, list_of_outputs @@ -306,7 +342,7 @@ def get_optical_spectrum(ensemble, w_array = None): Error, the configuration {} has no 'epsilon' data. """.format(i) raise ValueError(ERR) - + data = np.array(ensemble.all_properties[i]['epsilon']) if w_data is None: @@ -330,7 +366,7 @@ def get_optical_spectrum(ensemble, w_array = None): kind = 'cubic', bounds_error = False, fill_value = 'extrapolate') eps_imag = f_imag(w_array) w_data = w_array - + # Build the complex epsilon epsilon = eps_real + 1j * eps_imag @@ -338,5 +374,3 @@ def get_optical_spectrum(ensemble, w_array = None): n = np.sqrt(epsilon) return w_data, n - - diff --git a/Modules/Cluster.py b/Modules/Cluster.py index 3d12d2bd..a1efeef7 100644 --- a/Modules/Cluster.py +++ b/Modules/Cluster.py @@ -51,7 +51,7 @@ __CLUSTER_VNODES__ = "v_nodes" __CLUSTER_NNODES__ = "n_nodes" __CLUSTER_UNODES__ = "use_nodes" -__CLUSTER_VCPU__ = "v_cpu" +__CLUSTER_VCPU__ = "v_cpu" __CLUSTER_NCPU__ = "n_cpu" __CLUSTER_UCPU__ = "use_cpu" __CLUSTER_VTIME__ = "v_time" @@ -97,7 +97,7 @@ __CLUSTER_UPART__, __CLUSTER_VQOS__, __CLUSTER_NQOS__, __CLUSTER_UQOS__, __CLUSTER_INITSCRIPT__, __CLUSTER_MAXRECALC__, __CLUSTER_BATCHSIZE__, __CLUSTER_LOCALWD__, __CLUSTER_VACCOUNT__, __CLUSTER_UACCOUNT__, __CLUSTER_SSHCMD__, - __CLUSTER_SCPCMD__, __CLUSTER_WORKDIR__, __CLUSTER_TIMEOUT__, + __CLUSTER_SCPCMD__, __CLUSTER_WORKDIR__, __CLUSTER_TIMEOUT__, __CLUSTER_JOBNUMBER__, __CLUSTER_NPARALLEL__, __CLUSTER_NPOOLS__, __CLUSTER_ATTEMPTS__, __CLUSTER_PORT__] @@ -113,13 +113,13 @@ def parse_symbols(string): In a string that must be used on a shell command, replace all the special symbols in a way that they are correctly passed through the SHELL. - for example $USER => \$USER - + for example $USER => \$USER + Parameters ---------- string : str The string to be parsed - + Results ------- output : str @@ -134,7 +134,7 @@ def parse_symbols(string): class Cluster(object): - + def __init__(self, hostname=None, pwd=None, extra_options="", workdir = "", account_name = "", partition_name = "", qos_name = "", binary="pw.x -npool NPOOL -i PREFIX.pwi > PREFIX.pwo", mpi_cmd=r"srun --mpi=pmi2 -n NPROC"): @@ -143,17 +143,17 @@ def __init__(self, hostname=None, pwd=None, extra_options="", workdir = "", ================= It is strongly suggested to use public/private keys. However, sometimes for some reasons it is not possible to do so. - Then you must install sshpass, and provide the password + Then you must install sshpass, and provide the password for establishing the connection. - + Parameters ---------- hostname: The name of the host toward you want to connect. - E.G. + E.G. pippo@login.cineca.marconi.it pwd: - The password for the connection. + The password for the connection. If you use a private key, do not pass this argument. In the other case you must have sshpass package to allow the password communication. @@ -170,7 +170,7 @@ def __init__(self, hostname=None, pwd=None, extra_options="", workdir = "", qos_name: The QOS name for the job submission """ - + self.hostname = hostname self.pwd = pwd @@ -183,14 +183,14 @@ def __init__(self, hostname=None, pwd=None, extra_options="", workdir = "", res = os.system("sshpass > /dev/null") if res != 0: raise ValueError("Error, sshpass command not found, required to connect through password to server.") - + self.sshcmd = "sshpass -p '" + pwd + "' ssh " + extra_options self.scpcmd = "sshpass -p '" + pwd + "' scp " + extra_options.replace("-p", "-P") else: self.sshcmd = "ssh " + extra_options self.scpcmd = "scp " + extra_options.replace("-p", "-P") - + self.account_name = account_name self.partition_name = partition_name if partition_name: @@ -198,11 +198,11 @@ def __init__(self, hostname=None, pwd=None, extra_options="", workdir = "", self.qos_name = qos_name if qos_name: self.use_qos = True - + self.binary = binary self.mpi_cmd = mpi_cmd - self.workdir=r"" + self.workdir="" self.submit_command="sbatch" self.submit_name="SBATCH" self.terminal="/bin/bash" @@ -246,39 +246,39 @@ def __init__(self, hostname=None, pwd=None, extra_options="", workdir = "", # This is a set of lines to be runned before the calculation # It can be used to load the needed modules in the cluster self.load_modules="" - + # Setup the number of cpu, nodes and pools for the calculation self.n_nodes = 1 self.n_cpu = 1 self.n_pool = 1 - + # This is the default label, change it if you are going to submit # two different calculations in the same working directory self.label = "ESP_" - + # This is the maximum number of resubmissions after the expected one # from the batch size. It can be use to resubmit the failed jobs. self.max_recalc = 10 self.connection_attempts = 1 - + # This is the time string. Faster job will be the less in the queue, self.time="00:02:00" # 2minutes - + # The ram required for the calculation self.ram="10000Mb" # 10Gb - + # The default partition in which to submit calculations self.partition_name = "" - + # The default qos in which to submit calculations self.qos_name = "" # Still unused self.prefix_name = "prefix" # Variable in the calculator for differentiating the calculations - + # This directory is used to work with clusters. self.local_workdir = "cluster_work/" - + # The batch size is the maximum number of job to be submitted together. # The new jobs will be submitted only after a batch is compleated # Useful if the cluster has a limit for the maximum number of jobs allowed. @@ -309,7 +309,7 @@ def __setattr__(self, name, value): It will raise an exception if the attribute does not exists (with a suggestion of similar entries) """ - + if "fixed_attributes" in self.__dict__: if name in self.__total_attributes__: super(Cluster, self).__setattr__(name, value) @@ -331,17 +331,17 @@ def __setattr__(self, name, value): super(Cluster, self).__setattr__(name, value) - - - + + + def ExecuteCMD(self, cmd, raise_error = False, return_output = False, on_cluster = False): """ EXECUTE THE CMD ON THE CLUSTER ============================== - + This subroutine execute the cmd in the cluster, with the specified number of attempts. - + Parameters ---------- cmd: string @@ -349,11 +349,11 @@ def ExecuteCMD(self, cmd, raise_error = False, return_output = False, on_cluster raise_error : bool, optional If True (default) raises an error upon failure. return_output : bool, optional - If True (default False) the output of the command is + If True (default False) the output of the command is returned as second value. on_cluster : bool If true, the command is executed directly on the cluster through ssh - + Returns ------- success : bool @@ -366,7 +366,7 @@ def ExecuteCMD(self, cmd, raise_error = False, return_output = False, on_cluster if on_cluster: cmd = self.sshcmd + " {} '{}'".format(self.hostname, cmd) - + success = False output = "" for i in range(self.connection_attempts): @@ -381,7 +381,7 @@ def ExecuteCMD(self, cmd, raise_error = False, return_output = False, on_cluster sys.stderr.write(e) sys.stderr.write("\n") sys.stderr.flush() - + output = output.strip() status = p.wait() @@ -389,7 +389,7 @@ def ExecuteCMD(self, cmd, raise_error = False, return_output = False, on_cluster print('THREAD {} EXECUTE COMMAND: {}'.format(threading.get_native_id(), cmd)) if not err is None: sys.stderr.write(err) - + if status != 0: sys.stderr.write("Error with cmd: "+ cmd + "\n") sys.stderr.write(output + "\n") @@ -397,17 +397,17 @@ def ExecuteCMD(self, cmd, raise_error = False, return_output = False, on_cluster else: success = True break - + if raise_error and not success: raise IOError("Error while communicating with the cluster. More than %d attempts failed." % (i+1)) - - + + if return_output: return success, output return success - - - + + + def set_timeout(self, timeout): """ @@ -422,22 +422,22 @@ def set_timeout(self, timeout): self.use_timeout = True self.timeout = timeout - - - + + + def CheckCommunication(self): """ CHECK IF THE SERVER IS AVAILABLE ================================ - - This function return true if the server respond correctly, + + This function return true if the server respond correctly, false otherwise. """ - + cmd = self.sshcmd + " %s 'echo ciao'" % self.hostname - + status, output = self.ExecuteCMD(cmd, return_output = True) - + # print cmd # p = subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True) # output, err = p.communicate() @@ -449,14 +449,14 @@ def CheckCommunication(self): # sys.stderr.write("Error, cmd: " + cmd + "\n") # sys.stderr.write("Exit status:" + str(status)) # return False - + if output != "ciao": sys.stderr.write("Code exited correctly but echo did not answer properly\n") sys.stderr.write("cmd: " + cmd + "\n") sys.stderr.write("out: " + output + "\n") sys.stderr.write("expected " + "ciao" + "\n") return False - + return True def create_submission_script(self, labels): @@ -469,9 +469,9 @@ def create_submission_script(self, labels): Parameters ---------- - labels : list + labels : list It is a list of the labels of the calculations to be done. - + Returns ------- submission_header : string @@ -480,8 +480,8 @@ def create_submission_script(self, labels): # prepare the submission script submission = "#!" + self.terminal + "\n" - - + + # Add the submission options if self.use_nodes: @@ -511,14 +511,14 @@ def create_submission_script(self, labels): else: submission += "#{} {}={}\n".format(self.submit_name, adder_string, self.custom_params[add_parameter]) - + # Add the set -x option if self.add_set_minus_x: submission += "set -x\n" - + # Add the loading of the modules submission += self.load_modules + "\n" - + # Go to the working directory submission += "cd " + self.workdir + "\n" @@ -529,9 +529,9 @@ def create_submission_script(self, labels): other_input, other_output = self.additional_script_parameters(labels) submission += other_input - + # Use the xargs trick - #submission += "xargs -d " + r"'\n'" + " -L1 -P%d -a %s -- bash -c\n" % (n_togheder, + #submission += "xargs -d " + r"'\n'" + " -L1 -P%d -a %s -- bash -c\n" % (n_togheder, for i, lbl in enumerate(labels): submission += self.get_execution_command(lbl) @@ -560,7 +560,7 @@ def get_execution_command(self, label): # Get the MPI command replacing NPROC new_mpicmd = self.mpi_cmd.replace("NPROC", str(self.n_cpu)) - + # Replace the NPOOL variable and the PREFIX in the binary binary = self.binary.replace("NPOOL", str(self.n_pool)).replace("PREFIX", label) @@ -575,10 +575,10 @@ def prepare_input_file(self, structures, calc, labels): PREPARE THE INPUT FILE ====================== - This is specific for quantum espresso and must be inherit and replaced for + This is specific for quantum espresso and must be inherit and replaced for other calculators. - This crates the input files in the local working directory + This crates the input files in the local working directory self.local_workdir and it returns the list of all the files generated. @@ -606,7 +606,7 @@ def prepare_input_file(self, structures, calc, labels): for i, (label, structure) in enumerate(zip(labels, structures)): # Avoid thread conflict self.lock.acquire() - + try: calc.set_directory(self.local_workdir) calc.set_label(label) @@ -620,9 +620,9 @@ def prepare_input_file(self, structures, calc, labels): input_file = '{}.pwi'.format(label) output_file = '{}.pwo'.format(label) - list_of_inputs.append(input_file) + list_of_inputs.append(input_file) list_of_outputs.append(output_file) - except Exception as e: + except Exception as e: MSG = ''' Error while writing input file {}. @@ -633,9 +633,9 @@ def prepare_input_file(self, structures, calc, labels): # Release the lock on the threads - self.lock.release() - - + self.lock.release() + + return list_of_inputs, list_of_outputs def clean_localworkdir(self): @@ -684,12 +684,12 @@ def copy_files(self, list_of_input, list_of_output, to_server): for i, fname in enumerate(list_of_input): cmd += ' ' + fname - + # Remove all the previous output files on cluster rm_cmd = 'rm -f {}; '.format(os.path.join(self.workdir, fname)) - + os.system(cmd) if not os.path.exists(tar_file): @@ -706,7 +706,7 @@ def copy_files(self, list_of_input, list_of_output, to_server): if os.path.exists(os.path.join(self.local_workdir, fname)): os.remove(os.path.join(self.local_workdir, fname)) - + # Clean eventually input/output file of this very same calculation cmd = self.sshcmd + " %s '%s'" % (self.hostname, rm_cmd) self.ExecuteCMD(cmd, False) @@ -715,7 +715,7 @@ def copy_files(self, list_of_input, list_of_output, to_server): # print "Error while executing:", cmd # print "Return code:", cp_res # sys.stderr.write(cmd + ": exit with code " + str(cp_res) + "\n") -# +# # Copy the file into the cluster cmd = self.scpcmd + " %s %s:%s/" % (tar_file, self.hostname, self.workdir) @@ -746,7 +746,7 @@ def copy_files(self, list_of_input, list_of_output, to_server): tar_command = 'tar cf {} '.format(tar_name) for output in list_of_output: tar_command += ' ' + output - + compress_cmd = 'cd {}; {}'.format(self.workdir, tar_command) cmd = self.sshcmd + " %s '%s'" % (self.hostname, compress_cmd) @@ -755,7 +755,7 @@ def copy_files(self, list_of_input, list_of_output, to_server): print ("Error while compressing the outputs:", cmd, list_of_output, "\nReturn code:", cp_res) #return cp_res - + # Copy the tar and unpack cmd = self.scpcmd + "%s:%s %s/" % (self.hostname, os.path.join(self.workdir, tar_name), self.local_workdir) cp_res = self.ExecuteCMD(cmd, False) @@ -775,8 +775,8 @@ def copy_files(self, list_of_input, list_of_output, to_server): return cp_res - - + + def submit(self, script_location): """ @@ -789,33 +789,33 @@ def submit(self, script_location): return self.ExecuteCMD(cmd, True, return_output=True) - + Parameters ---------- script_localtion : string Path to the submission script inside the cluster. - + Results ------- success : bool - Result of the execution of the submission command. + Result of the execution of the submission command. It is what returned from self.ExecuteCMD(cmd, False) """ - - cmd = "{ssh} {host} '{submit_cmd} {script}'".format(ssh = self.sshcmd, host = self.hostname, - submit_cmd = self.submit_command, script = script_location) - if self.use_active_shell: - cmd = "{ssh} {host} -t '{shell} --login -c \"{submit_cmd} {script}\"'".format(ssh = self.sshcmd, - host = self.hostname, - submit_cmd = self.submit_command, script = script_location, + + cmd = "{ssh} {host} '{submit_cmd} {script}'".format(ssh = self.sshcmd, host = self.hostname, + submit_cmd = self.submit_command, script = script_location) + if self.use_active_shell: + cmd = "{ssh} {host} -t '{shell} --login -c \"{submit_cmd} {script}\"'".format(ssh = self.sshcmd, + host = self.hostname, + submit_cmd = self.submit_command, script = script_location, shell = self.terminal) - - - #cmd = self.sshcmd + " %s '%s %s/%s.sh'" % (self.hostname, self.submit_command, + + + #cmd = self.sshcmd + " %s '%s %s/%s.sh'" % (self.hostname, self.submit_command, # self.workdir, label+ "_" + str(indices[0])) - + return self.ExecuteCMD(cmd, False, return_output=True) def get_output_path(self, label): @@ -823,7 +823,7 @@ def get_output_path(self, label): Given the label of the submission, retrive the path of all the output files of that calculation """ - out_filename = os.path.join(self.workdir, label + ".pwo") + out_filename = os.path.join(self.workdir, label + ".pwo") return [out_filename] def read_results(self, calc, label): @@ -845,21 +845,21 @@ def read_results(self, calc, label): return results - def batch_submission(self, list_of_structures, calc, indices, + def batch_submission(self, list_of_structures, calc, indices, in_extension, out_extension, label = "ESP", n_togheder=1): """ BATCH SUBMISSION ================ - + This is a different kind of submission, it exploits xargs to perform a parallel submission of serveral structures in one single job. This is very good to avoid overloading the queue system manager, or when a limited number of jobs per user are allowed. - + NOTE: the number of structure in the list must be a divisor of the total number of processors. - + Parameters ---------- list_of_structures : list @@ -876,21 +876,21 @@ def batch_submission(self, list_of_structures, calc, indices, label : string, optional The root of the input file. n_togheder : int, optional (DO NOT USE) - If present, the job will lunch a new job immediately after the other - is ended. This is usefull to further reduce the number of submitted + If present, the job will lunch a new job immediately after the other + is ended. This is usefull to further reduce the number of submitted jobs. - + Results ------- list_of_results. Returns a list of results dicts, one for each structure. """ N_structs = len(list_of_structures) - + if n_togheder != 1: raise NotImplementedError("Error, n_togheder != 1 does not work!") - - + + # Prepare the input atoms app_list = "" new_ncpu = self.n_cpu * n_togheder @@ -898,21 +898,21 @@ def batch_submission(self, list_of_structures, calc, indices, submitted = [] submission_labels = [] - + for i in range(N_structs): # Prepare a typical label lbl = label + "_" + str(indices[i]) submission_labels.append(lbl) submitted.append(i) - - - # Create the input files + + + # Create the input files input_files, output_files = self.prepare_input_file(list_of_structures, calc, submission_labels) # Create the submission script submission = self.create_submission_script(submission_labels) - + # Add the submission script to the input files sub_name = label + "_" + str(indices[0]) + '.sh' sub_fpath = os.path.join(self.local_workdir, sub_name) @@ -920,16 +920,16 @@ def batch_submission(self, list_of_structures, calc, indices, f.write(submission) f.close() input_files.append(sub_name) - + if not self.copy_files(input_files, output_files, to_server = True): # Submission failed return submitted, indices, label - - + + # Run the simulation sub_script_loc = os.path.join(self.workdir, sub_name) cp_res, submission_output = self.submit(sub_script_loc) - + if self.nonblocking_command: job_id = self.get_job_id_from_submission_output(submission_output) @@ -948,11 +948,11 @@ def batch_submission(self, list_of_structures, calc, indices, else: print('[SUBMISSION {}] GOT OUTPUT'.format(threading.get_native_id())) - + return submitted, indices, label def collect_results(self, calc, submitted, indices, label): - """ + """ Collect back all the results from the submitted data. This operation needs to be performed in thread safe mode (all threads must be locked between this operation and when the results are actually assigned to the ensemble @@ -963,10 +963,10 @@ def collect_results(self, calc, submitted, indices, label): for i in submitted: # Prepare a typical label lbl = label + "_" + str(indices[i]) - + # Get the results try: - results[i] = self.read_results(calc, lbl) + results[i] = self.read_results(calc, lbl) except FileNotFoundError: sys.stderr.write("JOB {} | {} resulted in error:\n".format(i, lbl)) sys.stderr.write('File not found!\n') @@ -975,7 +975,7 @@ def collect_results(self, calc, submitted, indices, label): sys.stderr.write("JOB {} | {} resulted in error:\n".format(i, lbl)) print(e, file=sys.stderr) sys.stderr.flush() - + print('[SUBMISSION {}] GOT RESULTS:\n{}'.format(threading.get_native_id(), results)) return results @@ -984,7 +984,7 @@ def get_job_id_from_submission_output(self, output): """ GET THE JOB ID - Retreive the job id from the output of the submission. + Retreive the job id from the output of the submission. This depends on the software employed. It works for slurm. Returns None if the output contains an error @@ -995,12 +995,12 @@ def get_job_id_from_submission_output(self, output): try: id = output.split()[-1] - print('Understanding output "{}" => {}'.format(output, id)) + print('Understanding output "{}" => {}'.format(output, id)) return id except: print("Error, expected a standard output, but the result of the submission was: {}".format(output)) return None - + def check_job_finished(self, job_id, verbose = True): """ Check if the job identified by the job_id is finished @@ -1023,7 +1023,7 @@ def check_job_finished(self, job_id, verbose = True): sys.stderr.write("{}/{}/{} - {}:{}:{} | job {}: No response from the server \n".format(now.year, now.month, now.day, now.hour, now.minute, now.second, job_id)) sys.stderr.flush() return False - + data = l.split() if len(data) == 0: if verbose: @@ -1037,7 +1037,7 @@ def check_job_finished(self, job_id, verbose = True): sys.stderr.write("{}/{}/{} - {}:{}:{} | job {} still running\n".format(now.year, now.month, now.day, now.hour, now.minute, now.second, job_id)) sys.stderr.flush() return False - + # If I'm here it means I did not find the job, but the command returned at least 1 line (so it was correctly executed). if verbose: now = datetime.datetime.now() @@ -1050,28 +1050,28 @@ def check_job_finished(self, job_id, verbose = True): sys.stderr.flush() return False - - def run_atoms(self, ase_calc, ase_atoms, label="ESP", + + def run_atoms(self, ase_calc, ase_atoms, label="ESP", in_extension = ".pwi", out_extension=".pwo", n_nodes = 1, n_cpu=1, npool=1): """ RUN ATOMS ON THE CLUSTER ======================== - + This function runs the given atoms in the cluster, using the ase_calculator. - Note: the ase_calc must be a FileIOCalculator. + Note: the ase_calc must be a FileIOCalculator. For now it works with quantum espresso. """ - + # Setup the calculator with the atomic structure ase_atoms.set_calculator(ase_calc) - + # Write the input file ase.io.write("%s/%s%s" % (self.local_workdir, label, in_extension),ase_atoms,**ase_calc.parameters) - + # prepare the submission script submission = self.terminal + "\n" - + # Add the submission options if self.use_nodes: submission += "#%s %s%d\n" % (self.submit_name, self.v_nodes, n_nodes) @@ -1101,18 +1101,18 @@ def run_atoms(self, ase_calc, ase_atoms, label="ESP", # Add the loading of the modules submission += self.load_modules + "\n" - + # Go to the working directory submission += "cd " + self.workdir + "\n" - + # Get the real calculation command mpicmd = self.mpi_cmd.replace("NPROC", str(n_cpu)) binary = self.binary.replace("NPOOL", str(npool)).replace("PREFIX", label) - + submission += mpicmd + " " + binary + "\n" - + # First of all clean eventually input/output file of this very same calculation - cmd = self.sshcmd + " %s 'rm -f %s/%s%s %s/%s%s'" % (self.hostname, + cmd = self.sshcmd + " %s 'rm -f %s/%s%s %s/%s%s'" % (self.hostname, self.workdir, label, in_extension, self.workdir, label, out_extension) self.ExecuteCMD(cmd, False) @@ -1121,7 +1121,7 @@ def run_atoms(self, ase_calc, ase_atoms, label="ESP", # print "Error while executing:", cmd # print "Return code:", cp_res # sys.stderr.write(cmd + ": exit with code " + str(cp_res)) -# +# # Copy the input files into the target directory f = open("%s/%s.sh" % (self.local_workdir, label), "w") f.write(submission) @@ -1133,7 +1133,7 @@ def run_atoms(self, ase_calc, ase_atoms, label="ESP", # print "Error while executing:", cmd # print "Return code:", cp_res # sys.stderr.write(cmd + ": exit with code " + str(cp_res)) -# +# cmd = self.scpcmd + " %s/%s%s %s:%s" % (self.local_workdir, label, in_extension, self.hostname, self.workdir) cp_res = self.ExecuteCMD(cmd, False) #cp_res = os.system(cmd) @@ -1142,7 +1142,7 @@ def run_atoms(self, ase_calc, ase_atoms, label="ESP", #print "Return code:", cp_res #sys.stderr.write(cmd + ": exit with code " + str(cp_res)) return - + # Run the simulation cmd = self.sshcmd + " %s '%s %s/%s.sh'" % (self.hostname, self.submit_command, self.workdir, label) self.ExecuteCMD(cmd, False) @@ -1151,10 +1151,10 @@ def run_atoms(self, ase_calc, ase_atoms, label="ESP", # print "Error while executing:", cmd # print "Return code:", cp_res # sys.stderr.write(cmd + ": exit with code " + str(cp_res)) - + # Get the response - cmd = self.scpcmd + " %s:%s/%s%s %s/" % (self.hostname, self.workdir, label, out_extension, - self.local_workdir) + cmd = self.scpcmd + " %s:%s/%s%s %s/" % (self.hostname, self.workdir, label, out_extension, + self.local_workdir) cp_res = self.ExecuteCMD(cmd, False) #cp_res = os.system(cmd) if not cp_res: @@ -1162,39 +1162,39 @@ def run_atoms(self, ase_calc, ase_atoms, label="ESP", print ("Return code:", cp_res) sys.stderr.write(cmd + ": exit with code " + str(cp_res)) return - + # Get the results ase_calc.set_label("%s/%s" % (self.local_workdir, label)) ase_calc.read_results() - + return ase_calc.results - + def setup_from_namelist(self, namelist): """ SETUP THE CLUSTER WITH AN INPUTFILE =================================== - + This method setup the cluster using a custom input file. The inputfile must have the same shape of QuantumESPRESSO ones. The information about the cluster must be located in a namespace called as __CLUSTER_NAMELIST - + Parameters ---------- - namelist: + namelist: The parsed namelist dictionary. """ # Parse the namelist if needed if isinstance(namelist, str): namelist = CC.Methods.read_namelist(namelist) - + # Check if the cluster namelist is present if not __CLUSTER_NAMELIST__ in namelist.keys(): raise ValueError("Error, the parsed namelist must contain %s" % __CLUSTER_NAMELIST__) - + c_info = namelist[__CLUSTER_NAMELIST__] keys = c_info.keys() - + # Check if there is an unknown key for k in keys: if not k in __CLUSTER_KEYS__: @@ -1202,7 +1202,7 @@ def setup_from_namelist(self, namelist): s = "Did you mean something like:" + str( difflib.get_close_matches(k, __CLUSTER_KEYS__)) print (s) raise IOError("Error in cluster namespace: key '" + k +"' not recognized.\n" + s) - + # First of all, check if a template is present if __CLUSTER_TEMPLATE__ in keys: # Check if the environment variable has been defined @@ -1217,7 +1217,7 @@ def setup_from_namelist(self, namelist): if not os.path.exists(newfname): print ("Error, Environ variable", __TEMPLATE_ENV__, "exists, but no file", fname, "found.") raise IOError("Error while reading the cluster template.") - + self.setup_from_namelist(CC.Methods.read_namelist(newfname)) else: sys.stderr.write("Error, Environ variable" + __TEMPLATE_ENV__ + "exists, but not the directory it is pointing to.\n") @@ -1230,26 +1230,26 @@ def setup_from_namelist(self, namelist): for req_key in __CLUSTER_RKW__: if not req_key in keys: raise IOError("Error, the cluster configuration namelist requires the keyword: '" + req_key + "'") - - + + # Setup all the info if __CLUSTER_HOST__ in keys: self.hostname = c_info[__CLUSTER_HOST__] #print "HOST:", c_info[__CLUSTER_HOST__] - + if __CLUSTER_SSHCMD__ in keys: self.sshcmd = c_info[__CLUSTER_SSHCMD__] if __CLUSTER_SCPCMD__ in keys: self.scpcmd = c_info[__CLUSTER_SCPCMD__] - + if __CLUSTER_PWD__ in keys: self.pwd = c_info[__CLUSTER_PWD__] # Check if sshpass is present - res = os.system("sshpass > /dev/null") + res = os.system("sshpass > /dev/null") if res != 0: raise ValueError("Error, sshpass command not found, required to connect through password to server.") - + self.sshcmd = "sshpass -p '" + self.pwd + "' " + self.sshcmd self.scpcmd = "sshpass -p '" + self.pwd + "' " + self.scpcmd @@ -1258,13 +1258,13 @@ def setup_from_namelist(self, namelist): # Check if the password has been setup self.sshcmd += " -p {:.0f}".format(c_info[__CLUSTER_PORT__]) self.scpcmd += " -P {:.0f}".format(c_info[__CLUSTER_PORT__]) - + if __CLUSTER_ACCOUNT__ in keys: self.account_name = c_info[__CLUSTER_ACCOUNT__] - + if __CLUSTER_MPICMD__ in keys: self.mpi_cmd = c_info[__CLUSTER_MPICMD__] - + if __CLUSTER_TIMEOUT__ in keys: self.use_timeout = True self.timeout = int(c_info[__CLUSTER_TIMEOUT__]) @@ -1274,18 +1274,18 @@ def setup_from_namelist(self, namelist): if self.job_number < 1: print ("Error, the number of job per batch must be >= 1") raise ValueError("Error in the %s input variable." % __CLUSTER_JOBNUMBER__) - + if __CLUSTER_NPARALLEL__ in keys: self.n_together_def = int(c_info[__CLUSTER_NPARALLEL__]) - + if self.n_together_def < 1: print ("Error, the number of parallel jobs must be >= 1") raise ValueError("Error in the %s input variable." % __CLUSTER_NPARALLEL__) if self.n_together_def > self.job_number: print ("Error, the number of parallel runs must be <= than the number of runs per job") raise ValueError("Error, check the cluster keys %s and %s" % (__CLUSTER_NPARALLEL__, __CLUSTER_JOBNUMBER__)) - - + + if __CLUSTER_TERMINAL__ in keys: self.terminal = "#!" + c_info[__CLUSTER_TERMINAL__] if __CLUSTER_SUBCMD__ in keys: @@ -1337,10 +1337,10 @@ def setup_from_namelist(self, namelist): self.v_account = c_info[__CLUSTER_VACCOUNT__] if __CLUSTER_NPOOLS__ in keys: self.n_pool = int(c_info[__CLUSTER_NPOOLS__]) - + if __CLUSTER_ATTEMPTS__ in keys: self.connection_attempts = int(c_info[__CLUSTER_ATTEMPTS__]) - + if __CLUSTER_INITSCRIPT__ in keys: # Load the init script. # First lets parse the local environmental variables @@ -1348,84 +1348,84 @@ def setup_from_namelist(self, namelist): script_path = c_info[__CLUSTER_INITSCRIPT__] for key in k_env: script_path = script_path.replace("$" + key, os.environ[key]) - + script_path = os.path.abspath(script_path) - + # Check if the file exists if not os.path.exists(script_path): raise IOError("Error, the provided script path %s does not exists." % script_path) - + # Read the script and store the contnet f = open(script_path, "r") self.load_modules = f.read() f.close() - + if __CLUSTER_MAXRECALC__ in keys: self.max_recalc = int(c_info[__CLUSTER_MAXRECALC__]) if __CLUSTER_BATCHSIZE__ in keys: self.batch_size = int(c_info[__CLUSTER_BATCHSIZE__]) if __CLUSTER_LOCALWD__ in keys: wdir = c_info[__CLUSTER_LOCALWD__] - + # Parse the local environmental variables for ekey in os.environ.keys(): wdir = wdir.replace("$" + ekey, os.environ[ekey]) - + self.local_workdir = wdir - - + + if __CLUSTER_BINARY__ in keys: print ("Binary before parsing:") print (c_info[__CLUSTER_BINARY__]) self.binary = self.parse_string(c_info[__CLUSTER_BINARY__]) print ("Cluster binary setted to:") print (self.binary) - + # If all the cluster has been setted, setup the working directory if __CLUSTER_WORKDIR__ in keys: self.workdir = c_info[__CLUSTER_WORKDIR__] self.setup_workdir() - + def setup_workdir(self, verbose = True): """ SETUP THE WORKING DIRECTORY =========================== - + Parse the line contained in self.workdir in the claster to get working directory. It needs that the communication with the cluster has been correctly setted up. - + It will parse correctly environmental variables of the cluster. """ workdir = self.parse_string(self.workdir) - - sshcmd = self.sshcmd + " %s 'mkdir -p %s'" % (self.hostname, + + sshcmd = self.sshcmd + " %s 'mkdir -p %s'" % (self.hostname, workdir) - + self.ExecuteCMD(sshcmd, raise_error= True) -# +# # retval = os.system(sshcmd) # if retval != 0: # raise IOError("Error, while executing cmd: " + sshcmd) -# +# # if verbose: # print "Cluster workdir setted to:" # print workdir - + self.workdir = workdir - + def parse_string(self, string): """ PARSE STRING ============ - - Parse the given string on the cluster. + + Parse the given string on the cluster. It can be used to resolve environmental variables defined in the cluster. - + It will execute on the cluster the command: echo "string" and return the result of the cluster. - + Parameter --------- string : @@ -1436,63 +1436,63 @@ def parse_string(self, string): The same as input, but with the cluster environmental variables correctly parsed. """ - + # Open a pipe with the server # Use single ' to avoid string parsing by the local terminal cmd = "%s %s 'echo \"%s\"'" % (self.sshcmd, self.hostname, string) if self.use_active_shell: - cmd = "{ssh} {host} -t '{shell} --login -c \"echo {string}\"'".format(ssh = self.sshcmd, - host = self.hostname, - string = parse_symbols(string), + cmd = "{ssh} {host} -t '{shell} --login -c \"echo {string}\"'".format(ssh = self.sshcmd, + host = self.hostname, + string = parse_symbols(string), shell = self.terminal) #print cmd #print(cmd) - + status, output = self.ExecuteCMD(cmd, return_output = True, raise_error= True) -# +# # p = subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True) # output, err = p.communicate() # status = p.wait() # output = output.strip() # -# if not err is None: +# if not err is None: # sys.stderr.write(err) # if status != 0: # print "Command:", cmd # print "Error, status:", status # raise ValueError("Error, while connecting with the server.") - + return str(output) - + def compute_ensemble_batch(self, ensemble, cellconstructor_calc, get_stress = True, timeout=None): """ RUN THE ENSEMBLE WITH BATCH SUBMISSION ====================================== """ - + # Track the remaining configurations success = [False] * ensemble.N - + # Setup if the ensemble has the stress ensemble.has_stress = get_stress #ensemble.all_properties = [None] * ensemble.N - + # Check if the working directory exists if not os.path.isdir(self.local_workdir): os.makedirs(self.local_workdir) - - + + # Get the expected number of batch num_batch_offset = int(ensemble.N / self.batch_size) - + def compute_single_jobarray(jobs_id, calc): structures = [ensemble.structures[i].copy() for i in jobs_id] n_together = min(len(structures), self.n_together_def) subs, indices, labels = self.batch_submission(structures, calc, jobs_id, ".pwi", ".pwo", "ESP", n_together) - + # Thread safe operation self.lock.acquire() print("[THREAD {}] submitted calculations: {}".format(threading.get_native_id(), indices)) @@ -1504,7 +1504,7 @@ def compute_single_jobarray(jobs_id, calc): if res is None: continue - + # Check if the run was good check_e = "energy" in res check_f = "forces" in res @@ -1517,7 +1517,7 @@ def compute_single_jobarray(jobs_id, calc): print("ERROR IDENTIFYING STRUCTURE!") MSG = """ Error in thread {}. - Displacement between the expected structure {} + Displacement between the expected structure {} and the one readed from the calculator is of {} A. """.format(threading.get_native_id(), jobs_id[i], error_struct) @@ -1529,14 +1529,14 @@ def compute_single_jobarray(jobs_id, calc): continue else: print("[WARNING] no check on the structure.") - + is_success = check_e and check_f if get_stress: is_success = is_success and check_s - + if not is_success: continue - + res_only_extra = {x : res[x] for x in res if x not in ["energy", "forces", "stress", "structure"]} ensemble.all_properties[num].update(res_only_extra) ensemble.energies[num] = res["energy"] / units["Ry"] @@ -1557,7 +1557,7 @@ def compute_single_jobarray(jobs_id, calc): success[num] = is_success self.lock.release() - + # Run until some work has not finished recalc = 0 self.lock = threading.Lock() @@ -1566,17 +1566,17 @@ def compute_single_jobarray(jobs_id, calc): print("[CYCLE] SUCCESS: ", success) print("[CYCLE] STOPPING CONDITION:", np.sum(np.array(success, dtype = int) - 1)) - + # Get the remaining jobs false_mask = np.array(success) == False false_id = np.arange(ensemble.N)[false_mask] - + count = 0 # Submit in parallel jobs = [false_id[i : i + self.job_number] for i in range(0, len(false_id), self.job_number)] # Create a local copy of the calculator for each thread, to avoid conflicting modifications calculators = [cellconstructor_calc.copy() for i in range(0, len(jobs))] - + for k_th, job in enumerate(jobs): # Submit only the batch size if count >= self.batch_size: @@ -1585,11 +1585,11 @@ def compute_single_jobarray(jobs_id, calc): t.start() threads.append(t) count += 1 - + # Wait until all the job have finished for t in threads: t.join(timeout) - + print("[CYCLE] [END] SUCCESS: ", success) print("[CYCLE] [END] STOPPING CONDITION:", np.sum(np.array(success, dtype = int) - 1)) @@ -1598,43 +1598,43 @@ def compute_single_jobarray(jobs_id, calc): print ("Expected batch ordinary resubmissions:", num_batch_offset) raise ValueError("Error, resubmissions exceeded the maximum number of %d" % self.max_recalc) break - + print("CALCULATION ENDED: all properties: {}".format(ensemble.all_properties)) - - - + + + def compute_ensemble(self, ensemble, ase_calc, get_stress = True, timeout=None): """ RUN THE WHOLE ENSEMBLE ON THE CLUSTER ===================================== - + Parameters ---------- ensemble : The ensemble to be runned. """ - + # Check if the compute_ensemble batch must be done #if self.job_number != 1: self.compute_ensemble_batch(ensemble, ase_calc, get_stress, timeout) return - + """ # Track the remaining configurations success = [False] * ensemble.N - + # Setup if the ensemble has the stress ensemble.has_stress = get_stress - + # Check if the working directory exists if not os.path.isdir(self.local_workdir): os.makedirs(self.local_workdir) - + # Prepare the function for the simultaneous submission def compute_single(num, calc): atm = ensemble.structures[num].get_ase_atoms() - res = self.run_atoms(calc, atm, self.label + str(num), - n_nodes = self.n_nodes, + res = self.run_atoms(calc, atm, self.label + str(num), + n_nodes = self.n_nodes, n_cpu=self.n_cpu, npool = self.n_pool) if res: @@ -1654,19 +1654,19 @@ def compute_single(num, calc): # Remember, ase has a very strange definition of the stress ensemble.stresses[num, :, :] = -stress * units["Bohr"]**3 / units["Ry"] success[num] = True - + # Get the expected number of batch num_batch_offset = int(ensemble.N / self.batch_size) - + # Run until some work has not finished recalc = 0 while np.sum(np.array(success, dtype = int) - 1) != 0: threads = [] - + # Get the remaining jobs false_mask = np.array(success) == False false_id = np.arange(ensemble.N)[false_mask] - + count = 0 # Submit in parallel for i in false_id: @@ -1677,16 +1677,15 @@ def compute_single(num, calc): t.start() threads.append(t) count += 1 - - + + # Wait until all the job have finished for t in threads: t.join(timeout) - + recalc += 1 if recalc > num_batch_offset + self.max_recalc: print ("Expected batch ordinary resubmissions:", num_batch_offset) raise ValueError("Error, resubmissions exceeded the maximum number of %d" % self.max_recalc) break """ - diff --git a/Modules/Ensemble.py b/Modules/Ensemble.py index c06bbb9f..a14c55f4 100644 --- a/Modules/Ensemble.py +++ b/Modules/Ensemble.py @@ -4,7 +4,7 @@ import warnings import numpy as np import time - +#from scipy.special import tanh, sinh, cosh """ @@ -80,7 +80,7 @@ __JULIA_ERROR__ = "" try: import julia, julia.Main - julia.Main.include(os.path.join(os.path.dirname(__file__), + julia.Main.include(os.path.join(os.path.dirname(__file__), "fourier_gradient.jl")) __JULIA_EXT__ = True except: @@ -224,7 +224,7 @@ def __init__(self, dyn0, T0, supercell = None, **kwargs): You can omit the supercell keyword when initializing the ensemble if you are not sure on the value. """.format(nq, supercell) - + # Flag to use the fourier gradient # If true, avoid to compute the forces in real space at all. self.fourier_gradient = __JULIA_EXT__ @@ -257,7 +257,7 @@ def __init__(self, dyn0, T0, supercell = None, **kwargs): # Get the extra quantities self.all_properties = [] - # Initialize the q grid and lattice + # Initialize the q grid and lattice # For the fourier transform self.q_grid = np.array(self.dyn_0.q_tot) / CC.Units.A_TO_BOHR self.r_lat = np.zeros((nat_sc, 3), dtype = np.float64) @@ -445,7 +445,7 @@ def init(self): # Get the original u_disps self.u_disps_original = np.reshape( self.xats - np.tile(self.supercell_structure.coords, (self.N, 1,1)), - (self.N, 3 * nat_sc), + (self.N, 3 * nat_sc), order = "C" ) self.u_disps = self.u_disps_original.copy() @@ -455,7 +455,7 @@ def init(self): self.u_disps, self.q_grid, self.itau, - self.r_lat, + self.r_lat, ) self.u_disps_original_qspace = self.u_disps_qspace.copy() @@ -475,8 +475,8 @@ def init(self): self.sscha_energies[:] = -julia.Main.multiply_vector_vector_fourier( self.forces_qspace / CC.Units.A_TO_BOHR, self.u_disps_qspace * CC.Units.A_TO_BOHR - ) * 0.5 # The conversion is useless, keep it for clarity - + ) * 0.5 # The conversion is useless, keep it for clarity + self.sscha_forces = julia.Main.vector_q2r( self.sscha_forces_qspace, self.q_grid, @@ -488,7 +488,7 @@ def init(self): self.forces_qspace = None self.u_disps_qspace = None - + self.q_grid = np.array(self.dyn_0.q_tot) / CC.Units.A_TO_BOHR nat_sc = self.supercell_structure.N_atoms self.r_lat = np.zeros((nat_sc, 3), dtype = np.float64) @@ -709,7 +709,7 @@ def load(self, data_dir, population, N, verbose = False, load_displacements = Tr self.has_stress = True else: self.has_stress = False - + if timer: timer.execute_timed_function(self.init) else: @@ -732,7 +732,7 @@ def load(self, data_dir, population, N, verbose = False, load_displacements = Tr raise IOError(ERROR_MSG) - + def load_from_calculator_output(self, directory, out_ext = ".pwo", timer=None): """ @@ -836,7 +836,7 @@ def load_from_calculator_output(self, directory, out_ext = ".pwo", timer=None): timer.execute_timed_function(self.init) else: self.init() - + def save(self, data_dir, population, use_alat = False): """ @@ -988,7 +988,7 @@ def save_extxyz(self, filename, append_mode = True): self.save_enhanced_xyz(filename, append_mode = append_mode) raise ImportError("Error, this function requires ASE installed") - + ase_structs = [] for i, s in enumerate(self.structures): @@ -999,7 +999,7 @@ def save_extxyz(self, filename, append_mode = True): calculator = ase.calculators.singlepoint.SinglePointCalculator(struct, energy = energy, forces = forces, stress = CC.Methods.transform_voigt(stress)) - + struct.set_calculator(calculator) ase_structs.append(struct) @@ -1220,8 +1220,8 @@ def load_bin(self, data_dir, population_id = 1, avoid_loading_dyn = False, timer timer.execute_timed_function(self.init) else: self.init() - - + + def init_from_structures(self, structures): """ Initialize the ensemble from the given list of structures @@ -1250,7 +1250,7 @@ def init_from_structures(self, structures): for i, s in enumerate(self.structures): # Get the displacements self.xats[i, :, :] = s.coords - + # Initialize the supercell super_struct, itau = self.dyn_0.structure.generate_supercell(self.supercell, get_itau=True) self.supercell_structure = super_struct @@ -1267,13 +1267,13 @@ def init_from_structures(self, structures): self.u_disps, self.q_grid, self.itau, - self.r_lat + self.r_lat ) self.u_disps_original_qspace = julia.Main.vector_r2q( self.u_disps_original, self.q_grid, self.itau, - self.r_lat + self.r_lat ) # Get the dynamical matrix in the fourier format nat = self.dyn_0.structure.N_atoms @@ -1291,7 +1291,7 @@ def init_from_structures(self, structures): self.sscha_energies[:] = julia.Main.multiply_vector_vector_fourier( self.sscha_forces_qspace, self.u_disps_original_qspace * CC.Units.A_TO_BOHR - ) * 0.5 + ) * 0.5 self.rho = np.ones(self.N, dtype = np.float64) @@ -1310,7 +1310,7 @@ def init_from_structures(self, structures): self.init_q_opposite() # Initialize the supercell - + self.q_grid = np.array(self.dyn_0.q_tot) / CC.Units.A_TO_BOHR nat_sc = self.supercell_structure.N_atoms self.r_lat = np.zeros((nat_sc, 3), dtype = np.float64) @@ -1321,13 +1321,13 @@ def init_from_structures(self, structures): self.init_q_opposite() - + def init_q_opposite(self): """ Setup the inverse q points. - + This subroutine identifies the inverse of each q point from the dynamical matrix""" if self.fourier_gradient: bg = self.current_dyn.structure.get_reciprocal_vectors() / (2* np.pi ) @@ -1644,23 +1644,23 @@ def update_weights_fourier(self, new_dynamical_matrix, newT, timer=None): super_structure = self.supercell_structure.copy() Nat_sc = super_structure.N_atoms - + #self.u_disps[:,:] = self.xats.reshape((self.N, 3*Nat_sc)) - np.tile(super_structure.coords.ravel(), (self.N,1)) #old_disps[:,:] = self.xats.reshape((self.N, 3*Nat_sc)) - np.tile(super_struct0.coords.ravel(), (self.N,1)) - - + + # Get the new displacements. # In fourier space this only affects the q=0 point self.update_displacements(new_dynamical_matrix.structure) # Get the lattice vectors nat_sc = super_structure.N_atoms - + # Get the displacements according to Fourier u_disp_fourier_new = self.u_disps_qspace u_disp_fourier_old = self.u_disps_original_qspace - + if changed_dyn: if timer: w_new, pols, wqn, polsqn = timer.execute_timed_function(new_dynamical_matrix.DiagonalizeSupercell, return_qmodes=True) @@ -1690,7 +1690,7 @@ def update_weights_fourier(self, new_dynamical_matrix, newT, timer=None): t2 = time.time() if timer: timer.add_timer("Prepare dynq for julia", t2-t1) - + # Get the forces [Ry/A] and energies [Ry] self.sscha_forces_qspace = - julia.Main.multiply_matrix_vector_fourier( dynq, @@ -1700,7 +1700,7 @@ def update_weights_fourier(self, new_dynamical_matrix, newT, timer=None): self.sscha_energies[:] = -julia.Main.multiply_vector_vector_fourier( self.sscha_forces_qspace / CC.Units.A_TO_BOHR, u_disp_fourier_new * CC.Units.A_TO_BOHR - ) * 0.5 + ) * 0.5 # Go back to real space and convert in Ry/Angstrom # self.sscha_forces = julia.Main.vector_q2r( @@ -1708,7 +1708,7 @@ def update_weights_fourier(self, new_dynamical_matrix, newT, timer=None): # np.array(new_dynamical_matrix.q_tot), # self.itau, # r_lat - # ) * CC.Units.A_TO_BOHR + # ) * CC.Units.A_TO_BOHR t2 = time.time() @@ -1732,7 +1732,7 @@ def update_weights_fourier(self, new_dynamical_matrix, newT, timer=None): w = np.array(w/2, dtype = np.float64) # Get the a_0 - old_a = SCHAModules.thermodynamic.w_to_a(w, self.T0) + old_a = self.w_to_a(w, self.T0) # Now do the same for the new dynamical matrix #new_super_dyn = new_dynamical_matrix.GenerateSupercellDyn(self.supercell) @@ -1770,7 +1770,7 @@ def update_weights_fourier(self, new_dynamical_matrix, newT, timer=None): w= w_new[~trans_mask] w = np.array(w/2, dtype = np.float64) - new_a = SCHAModules.thermodynamic.w_to_a(w, newT) + new_a = self.w_to_a(w, newT) # Get the new displacements in the supercell @@ -1814,7 +1814,7 @@ def update_weights_fourier(self, new_dynamical_matrix, newT, timer=None): self.pols_q_0, self.dyn_0.structure.get_masses_array(), np.float64(self.T0) - ) + ) t2 = time.time() if timer: @@ -1843,7 +1843,7 @@ def update_weights_fourier(self, new_dynamical_matrix, newT, timer=None): # Get the normalization ratio #norm = np.sqrt(np.abs(np.linalg.det(ups_new) / np.linalg.det(ups_old))) t2 = time.time() - self.rho = np.prod( old_a / new_a) * np.exp(-0.5 * (uYu_new - uYu_old) ) + self.rho = np.prod( old_a / new_a) * np.exp(-0.5 * (uYu_new - uYu_old) ) t3 = time.time() if timer: @@ -1862,7 +1862,7 @@ def update_weights_fourier(self, new_dynamical_matrix, newT, timer=None): # Get the new displacement #self.u_disps[i, :] = self.structures[i].get_displacement(new_super_dyn.structure).reshape(3 * new_super_dyn.structure.N_atoms) #self.u_disps[i, :] = (self.xats[i, :, :] - super_structure.coords).reshape( 3*Nat_sc ) - + if timer: self.current_dyn = timer.execute_timed_function(new_dynamical_matrix.Copy) #print( "Time elapsed to update weights the sscha energies, forces and displacements:", t1 - t3, "s") @@ -1908,7 +1908,7 @@ def update_weights(self, new_dynamical_matrix, newT, update_q = False, timer=Non self.u_disps[:,:] = self.xats.reshape((self.N, 3*Nat_sc)) - np.tile(super_structure.coords.ravel(), (self.N,1)) old_disps[:,:] = self.xats.reshape((self.N, 3*Nat_sc)) - np.tile(super_struct0.coords.ravel(), (self.N,1)) - + if changed_dyn: if timer: @@ -1954,7 +1954,7 @@ def update_weights(self, new_dynamical_matrix, newT, update_q = False, timer=Non w = np.array(w/2, dtype = np.float64) # Get the a_0 - old_a = SCHAModules.thermodynamic.w_to_a(w, self.T0) + old_a = self.w_to_a(w, self.T0) # Now do the same for the new dynamical matrix #new_super_dyn = new_dynamical_matrix.GenerateSupercellDyn(self.supercell) @@ -1992,7 +1992,7 @@ def update_weights(self, new_dynamical_matrix, newT, update_q = False, timer=Non w= w_new[~trans_mask] w = np.array(w/2, dtype = np.float64) - new_a = SCHAModules.thermodynamic.w_to_a(w, newT) + new_a = self.w_to_a(w, newT) # Get the new displacements in the supercell @@ -2178,14 +2178,14 @@ def get_fourier_forces(self, get_error=True): transform to avoid doing the average once more. """ - + delta_forces = np.real(self.forces_qspace[:, :, 0] - self.sscha_forces_qspace[:, :, 0]) nq = self.q_grid.shape[0] delta_forces /= np.sqrt(nq) sum_f = self.rho.dot(delta_forces) N_eff = np.sum(self.rho) f_average = sum_f / N_eff - + if get_error: sum_f2 = self.rho.dot(delta_forces**2) error_f = np.sqrt(sum_f2 / N_eff - f_average**2) / np.sqrt(N_eff) @@ -2453,19 +2453,19 @@ def get_fourier_gradient(self, timer=None): GET GRADIENT IN FOURIER SPACE ============================= - This subroutine computes the gradient in foureir space + This subroutine computes the gradient in fourier space - It employes the acceleration available thanks + It employes the acceleration available thanks to the julia code. In brief, the calculation performs: .. math:: - \left< u(q) \delta f(-q)\right> + \left< u(q) \delta f(-q)\right> - To get the gradient in the q space. - The method can be easily parallelized, + To get the gradient in the q space. + The method can be easily parallelized, since the julia subroutine returns also the average of the square. @@ -2473,7 +2473,7 @@ def get_fourier_gradient(self, timer=None): ---------- - timer : Timer object, optional If present, the time spent in the julia code is added to the timer. - + Results ------- - grad : ndarray (3*nat x 3*nat) @@ -2487,11 +2487,11 @@ def get_fourier_gradient(self, timer=None): Error while loading the julia module. This subroutine requires the julia speedup. install julia as specified in the guide, - and execute the python script using the python-jl - interpreter. + and execute the python script using the python-jl + interpreter. """ raise ImportError(MSG) - + t1 = time.time() # Check if the opposite q are initialize, otherwise do it once for all if self.q_opposite_index is None: @@ -2501,14 +2501,14 @@ def get_fourier_gradient(self, timer=None): nq = len(self.current_dyn.q_tot) nat = self.current_dyn.structure.N_atoms nat_sc = self.structures[0].N_atoms - + Y_qspace = julia.Main.get_upsilon_fourier( self.w_q_current, self.pols_q_current, self.current_dyn.structure.get_masses_array(), np.float64(self.current_T) ) - + t3 = time.time() # Get the v tilde (Yu) vector in fourier space (Bohr^-1) @@ -2522,15 +2522,15 @@ def get_fourier_gradient(self, timer=None): t4 = time.time() - phi_grad = np.zeros((nq, 3*nat, 3*nat), dtype = np.complex128) + phi_grad = np.zeros((nq, 3*nat, 3*nat), dtype = np.complex128) phi_grad2 = np.zeros((nq, 3*nat, 3*nat), dtype = np.float64) - + # if self.itau is not None: # super_struct = self.supercell_structure # itau = self.itau # else: # # Create the lattices - # if timer is not None: + # if timer is not None: # super_struct, itau = timer.execute_timed_function(self.current_dyn.structure.generate_supercell, self.supercell, get_itau=True) # else: # super_struct, itau = self.current_dyn.structure.generate_supercell(self.supercell, get_itau=True) @@ -2546,14 +2546,14 @@ def get_fourier_gradient(self, timer=None): # r_lat *= CC.Units.A_TO_BOHR # q_tot = np.array(self.current_dyn.q_tot) / CC.Units.A_TO_BOHR - + t5 = time.time() # f_vector = (self.forces - self.sscha_forces).reshape( (self.N, 3 * nat_sc)) * CC.Units.BOHR_TO_ANGSTROM # u_vector = self.u_disps.reshape( (self.N, 3 * nat_sc)) / CC.Units.BOHR_TO_ANGSTROM # Call the julia subroutine - delta_forces = self.forces_qspace - self.sscha_forces_qspace + delta_forces = self.forces_qspace - self.sscha_forces_qspace delta_forces /= CC.Units.A_TO_BOHR # Compute the gradient splitting the configurations @@ -2576,7 +2576,7 @@ def get_gradient_worker(cfgs, timer=None): result = np.concatenate((phi_grad, phi_grad2), axis=0) return result - + # Get the range of the configurations to be computed for each processor configs_ranges = CC.Settings.split_configurations(self.N) @@ -2592,7 +2592,7 @@ def get_gradient_worker(cfgs, timer=None): phi_grad = gradient_and_error[:nq, :, :] phi_grad2 = gradient_and_error[nq:, :, :] - + # Divide by the total weights n_tot = np.sum(self.rho) phi_grad /= n_tot @@ -2628,12 +2628,12 @@ def work_function(argument, timer=None): gradient, _ = new_ensemble.get_preconditioned_gradient(*args, timer=timer, **kwargs) av_ensemble = np.sum(new_ensemble.rho) - + return gradient * av_ensemble # Get the range of the configurations to be computed for each processor configs_ranges = CC.Settings.split_configurations(self.N) - + if timer: gradient = timer.execute_timed_function(CC.Settings.GoParallel, work_function, configs_ranges, "+") else: @@ -2643,10 +2643,10 @@ def work_function(argument, timer=None): return gradient, np.zeros_like(gradient) + 1 - - - + + + def get_preconditioned_gradient(self, subtract_sscha = True, return_error = False, use_ups_supercell = True, preconditioned = 1, @@ -2698,7 +2698,7 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals super_struct = self.current_dyn.structure.generate_supercell(self.supercell) #supercell_dyn = self.current_dyn.GenerateSupercellDyn(self.supercell) - + # Dyagonalize if timer: w, pols = timer.execute_timed_function(self.current_dyn.DiagonalizeSupercell) @@ -2712,7 +2712,7 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals ityp = super_struct.get_ityp() + 1 # Py to fortran convertion mass = np.array(list(super_struct.masses.values())) - + log_err = "err_yesrho" mass *= 2 @@ -2742,7 +2742,7 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals if fast_grad or not preconditioned: if timer: - grad, grad_err = timer.execute_timed_function(SCHAModules.get_gradient_supercell, + grad, grad_err = timer.execute_timed_function(SCHAModules.get_gradient_supercell, self.rho, u_disp, eforces, w, pols, trans, self.current_T, mass, ityp, log_err, self.N, nat, 3*nat, len(mass), preconditioned, @@ -2771,7 +2771,7 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals # Perform the fourier transform if return_error: # Check if a multiprocessing function can be exploited - if hasattr(CC.Phonons, 'GetDynQFromFCSupercell_parallel') and CC.Settings.GetNProc() > 1: + if hasattr(CC.Phonons, 'GetDynQFromFCSupercell_parallel') and CC.Settings.GetNProc() > 1: if timer: q_grad,q_grad_err = timer.execute_timed_function(CC.Phonons.GetDynQFromFCSupercell_parallel, grad, np.array(self.current_dyn.q_tot), @@ -2798,7 +2798,7 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals self.current_dyn.structure, super_struct) else: if timer: - q_grad = timer.execute_timed_function(CC.Phonons.GetDynQFromFCSupercell, + q_grad = timer.execute_timed_function(CC.Phonons.GetDynQFromFCSupercell, grad, np.array(self.current_dyn.q_tot), self.current_dyn.structure, super_struct) else: @@ -3654,7 +3654,7 @@ def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, v else: w, pols = w_pols - a = SCHAModules.thermodynamic.w_to_a(w, self.current_T) + a = self.w_to_a(w, self.current_T) n_modes = len(w) @@ -3716,7 +3716,7 @@ def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, v else: qe_sym.SetupFromSPGLIB() qe_sym.ApplySymmetryToTensor3(d3) - + if verbose: print("Saving the third order force constants as d3_realspace_sym.npy [Ha units]") @@ -3935,7 +3935,7 @@ def split(self, split_mask): N = np.sum(split_mask.astype(int)) ens = Ensemble(self.dyn_0, self.T0, self.dyn_0.GetSupercell()) ens.init_from_structures(structs) - + ens.force_computed[:] = self.force_computed[split_mask] ens.stress_computed[:] = self.stress_computed[split_mask] @@ -4195,8 +4195,16 @@ def get_energy_forces(self, ase_calculator, compute_stress = True, stress_numeri timer.execute_timed_function(self.init) else: self.init() + def w_to_a(self,w, T): + n = len(w) + a = np.zeros(n) + if T == 0.0: + a[:] = np.sqrt(1.0 / (2.0 * w)) + else: + a[:] = np.sqrt((1.0 / np.tanh(0.5 * w * 315774.65221921849 / T)) / (2.0 * w)) + return a - +#------------------------------------------------------------------------------- def _wrapper_julia_get_upsilon_q(*args, **kwargs): """Worker function, just for testing""" return julia.Main.get_upsilon_fourier(*args, **kwargs) @@ -4210,7 +4218,7 @@ def _wrapper_julia_vector_q2r(*args, **kwargs): Parameters ---------- - - vector_q : ndarray (3*nat, nq, n_random) + - vector_q : ndarray (3*nat, nq, n_random) - q_tot : ndarray (nq, 3) - itau : ndarray (nat_sc) - R_lat : ndarray (nat_sc, 3) @@ -4219,7 +4227,7 @@ def _wrapper_julia_vector_q2r(*args, **kwargs): ------- - vector_sc : ndarray (n_random, 3*nat_sc) - + """ return julia.Main.vector_q2r(*args, **kwargs) @@ -4234,14 +4242,14 @@ def _wrapper_julia_vector_r2q(*args, **kwargs): Parameters ---------- - - vector_sc : ndarray (n_random, 3*nat_sc) + - vector_sc : ndarray (n_random, 3*nat_sc) - q_tot : ndarray (nq, 3) - itau : ndarray (nat_sc) - R_lat : ndarray (nat_sc, 3) Returns ------- - + - vector_q : ndarray (3*nat, nq, n_random) """ @@ -4263,12 +4271,12 @@ def _wrapper_julia_matrix_vector_fourier(*args, **kwargs): Returns ------- - + - vector_q : ndarray (n_random, 3*nat, nq) The result of the matrix vector multiplication """ - return julia.Main.multiply_matrix_vector_fourier(*args, + return julia.Main.multiply_matrix_vector_fourier(*args, **kwargs) def _wrapper_julia_vector_vector_fourier(*args, **kwargs): @@ -4291,7 +4299,5 @@ def _wrapper_julia_vector_vector_fourier(*args, **kwargs): The result of the vector vector multiplication """ - return julia.Main.multiply_vector_vector_fourier(*args, + return julia.Main.multiply_vector_vector_fourier(*args, **kwargs) - - diff --git a/Modules/SchaMinimizer.py b/Modules/SchaMinimizer.py index 7be8844f..fe00975f 100644 --- a/Modules/SchaMinimizer.py +++ b/Modules/SchaMinimizer.py @@ -111,7 +111,7 @@ class SSCHA_Minimizer(object): def __init__(self, ensemble = None, root_representation = "normal", kong_liu_ratio = 0.5, meaningful_factor = 0.2, - minimization_algorithm = "sdes", lambda_a = 1, + minimization_algorithm = "sdes", lambda_a = 1, timer=None, **kwargs): """ This class create a minimizer to perform the sscha minimization. @@ -437,7 +437,7 @@ def minimization_step(self, custom_function_gradient = None, timer=None): timer.add_timer("Symmetrize in the supercell", t_7 - t_6, timer=timer_apply) timer.add_timer("Return in fourier space", t_8 - t_7, timer=timer_return) - + # Apply the sum rule at gamma if timer is not None: timer.execute_timed_function(CC.symmetries.CustomASR, dyn_grad[0,:,:]) @@ -614,7 +614,7 @@ def minimization_step(self, custom_function_gradient = None, timer=None): timer.execute_timed_function(self.minimizer.update_dyn, new_kl_ratio, dyn_grad, struct_grad) new_dyn, new_struct = timer.execute_timed_function(self.minimizer.get_dyn_struct) else: - self.minimizer.update_dyn(new_kl_ratio, dyn_grad, struct_grad) + self.minimizer.update_dyn(new_kl_ratio, dyn_grad, struct_grad) new_dyn, new_struct = self.minimizer.get_dyn_struct() @@ -1066,7 +1066,7 @@ def init(self, verbosity = False, delete_previous_data = True, init_timer = True If true, it will clean previous minimizations from the free energies, gradients... init_timer : bool If true, it will initialize the timer for timing the function. - This can slightly slowdown the minimization, + This can slightly slowdown the minimization, however enables the use of the timer to check the time spent in each step. """ @@ -1295,7 +1295,7 @@ def run(self, verbose = 1, custom_function_pre = None, custom_function_post = No # Perform the minimization step if timer is not None: - timer.execute_timed_function(self.minimization_step, custom_function_gradient) + timer.execute_timed_function(self.minimization_step, custom_function_gradient) im_freqs = timer.execute_timed_function(self.check_imaginary_frequencies) else: self.minimization_step(custom_function_gradient) diff --git a/requirements.txt b/requirements.txt index e0237666..b968921c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,5 +2,5 @@ setuptools numpy scipy ase -spglib +spglib<=2.2 julia diff --git a/requirements2.txt b/requirements2.txt index b77539ae..590aa382 100644 --- a/requirements2.txt +++ b/requirements2.txt @@ -2,4 +2,4 @@ setuptools numpy scipy ase==3.16.0 -spglib +spglib<=2.2 diff --git a/tests/aiida_ensemble/test_aiida_ensemble.py b/tests/aiida_ensemble/test_aiida_ensemble.py index 46125ed7..1ea70520 100644 --- a/tests/aiida_ensemble/test_aiida_ensemble.py +++ b/tests/aiida_ensemble/test_aiida_ensemble.py @@ -27,7 +27,7 @@ def test_clean_runs(): ensemble.force_computed = np.array([True, False, True, True], dtype=bool) ensemble.stress_computed = np.copy(ensemble.force_computed) ensemble._clean_runs() - + assert all(ensemble.force_computed) assert len(ensemble.force_computed) == 3 assert len(ensemble.stress_computed) == 3 @@ -44,25 +44,25 @@ def test_get_running_workchains(generate_workchain_pw_node): running = ProcessState.RUNNING workchains = [ - generate_workchain_pw_node(process_state=finished, exit_status=0, label='T_300_id_0'), + generate_workchain_pw_node(process_state=finished, exit_status=0, label='T_300_id_0'), generate_workchain_pw_node(process_state=finished, exit_status=300, label='T_300_id_1'), ] success = [False, False] - + wcs_left = get_running_workchains(workchains=workchains, success=success) - + assert not wcs_left assert success == [True, False] workchains = [ - generate_workchain_pw_node(process_state=running, label='T_300_id_0'), + generate_workchain_pw_node(process_state=running, label='T_300_id_0'), generate_workchain_pw_node(process_state=finished, label='T_300_id_1'), generate_workchain_pw_node(process_state=running, label='T_300_id_2'), ] success = [False, False, False] - + wcs_left = get_running_workchains(workchains=workchains, success=success) - + assert wcs_left == [workchains[0], workchains[2]] assert success == [False, True, False] @@ -72,15 +72,14 @@ def test_submit_and_get_workchains(fixture_code): """Test the :func:`sscha.aiida_ensemble.submit_and_get_workchains` method.""" from cellconstructor.Structure import Structure from sscha.aiida_ensemble import submit_and_get_workchains - + pw_code = fixture_code('quantumespresso.pw') structures = [Structure(nat=1) for _ in range(5)] - + workchains = submit_and_get_workchains( - structures=structures, + structures=structures, pw_code=pw_code, temperature=300, ) - + assert len(workchains) == 5 -