diff --git a/Modules/Classify.py b/Modules/Classify.py index f7b117c0..939744b3 100644 --- a/Modules/Classify.py +++ b/Modules/Classify.py @@ -436,10 +436,10 @@ def recognize_triplet(dyn, orbit1a, map_uc, symprec=1e-5, verbose=False): Maps an arbitrary triplet to its reference triplet, permutaion, and symmtry operation after the classification. Dimension [Natom,Natom_sc,Natom_sc,3] """ - start_time = time.time() + # start_time = time.time() - if verbose: - print("===== STARTING TRIPLET CLASSIFICATION =====") + # if verbose: + # print("===== STARTING TRIPLET CLASSIFICATION =====") # Get the symmetries of the unit cell spg_syms = spglib.get_symmetry(dyn.structure.get_spglib_cell(), symprec) @@ -536,21 +536,21 @@ def recognize_triplet(dyn, orbit1a, map_uc, symprec=1e-5, verbose=False): nref3 = SCHAModules.module_hess.get_nref3(nat,ntot,tot3,nsym,orbit1a,map_uc) orbit3a, orbit3s, norbit, indep_fc, n_indep_fc, kernel, mapping_triplet = SCHAModules.module_hess.recognize_triplet(nat,ntot,tot3,nref3,nsym,orbit1a,map_uc,nontrivial,M,verbose) - end_time = time.time() - execution_time = end_time - start_time + # end_time = time.time() + # execution_time = end_time - start_time - # Print detailed info - if verbose: - print(" ") - print("Total triplets:", tot3) - print("Number of Orbits:", nref3) - print("Total number of independent elements:", sum(n_indep_fc[:nref3])) - print(" ") + # # Print detailed info + # if verbose: + # print(" ") + # print("Total triplets:", tot3) + # print("Number of Orbits:", nref3) + # print("Total number of independent elements:", sum(n_indep_fc[:nref3])) + # print(" ") - print("execution_time in triplet recognition:", execution_time, " s") - print("===== TRIPLET CLASSIFICATION and GJ ELIMINATION FINISHED ======") - print(" ") + # print("execution_time in triplet recognition:", execution_time, " s") + # print("===== TRIPLET CLASSIFICATION and GJ ELIMINATION FINISHED ======") + # print(" ") sys.stdout.flush() @@ -603,10 +603,10 @@ def recognize_quadruplet(dyn, orbit1a, map_uc, verbose=False, symprec=1e-5): the classification. Dimension [Natom,Natom_sc,Natom_sc,3] """ - start_time = time.time() + # start_time = time.time() - if verbose: - print("===== STARTING QUADRUPLET CLASSIFICATION =====") + # if verbose: + # print("===== STARTING QUADRUPLET CLASSIFICATION =====") # Get the symmetries of the unit cell spg_syms = spglib.get_symmetry(dyn.structure.get_spglib_cell(), symprec) @@ -671,19 +671,19 @@ def recognize_quadruplet(dyn, orbit1a, map_uc, verbose=False, symprec=1e-5): nref4 = SCHAModules.module_hess.get_nref4(nat,ntot,tot4,nsym,orbit1a,map_uc) orbit4a, orbit4s, norbit, indep_fc, n_indep_fc, kernel, mapping_quadruplet = SCHAModules.module_hess.recognize_quadruplet(nat,ntot,tot4,nref4,nsym,orbit1a,map_uc,nontrivial,M,verbose) - end_time = time.time() - execution_time = end_time - start_time + # end_time = time.time() + # execution_time = end_time - start_time - if verbose: - print(" ") - print("Total quadruplets:", tot4) - print("Number of Orbits:",nref4) - print("Total number of independent elements:", sum(n_indep_fc[:nref4])) - print(" ") + # if verbose: + # print(" ") + # print("Total quadruplets:", tot4) + # print("Number of Orbits:",nref4) + # print("Total number of independent elements:", sum(n_indep_fc[:nref4])) + # print(" ") - print("execution_time in quadruplet recognition:", execution_time, " s") - print("===== QUADRUPLET CLASSIFICATION and GJ ELIMINATION FINISHED ======") - print(" ") + # print("execution_time in quadruplet recognition:", execution_time, " s") + # print("===== QUADRUPLET CLASSIFICATION and GJ ELIMINATION FINISHED ======") + # print(" ") return(orbit4a[:nref4], orbit4s[:nref4], norbit[:nref4], indep_fc[:nref4], n_indep_fc[:nref4], kernel[:nref4], Rot, mapping_quadruplet) # Small functions required in main subroutines. diff --git a/Modules/Ensemble.py b/Modules/Ensemble.py index 78f6d8af..8c1945cd 100644 --- a/Modules/Ensemble.py +++ b/Modules/Ensemble.py @@ -3869,7 +3869,7 @@ def get_free_energy_hessian_dev(self, include_v4 = False, do_scf = True, eps = 1 nat_sc = dyn_supercell.structure.N_atoms n_modes = nat_sc*3 - mapping, rot_cart, map_uc, map_tr, T_list, T_list_frac = Classify.map_singlet(self.current_dyn, verbose = verbose) + mapping, rot_cart, map_uc, map_tr, T_list, T_list_frac = Classify.map_singlet(self.current_dyn, verbose = False) orbit2a, orbit2s, norbit, indep_elem, n_indep_elem, tensor = Classify.recognize_doublet(self.current_dyn, mapping, map_uc, verbose = verbose) @@ -4712,5 +4712,4 @@ def _wrapper_julia_vector_vector_fourier(*args, **kwargs): """ return julia.Main.multiply_vector_vector_fourier(*args, - **kwargs) - \ No newline at end of file + **kwargs) \ No newline at end of file diff --git a/Modules/qClassify.py b/Modules/qClassify.py index deec18c4..ed98db19 100644 --- a/Modules/qClassify.py +++ b/Modules/qClassify.py @@ -51,25 +51,31 @@ def map_singlet(q_list, q_list_frac, rcell, rot_cart): q_list_frac_fixed = np.empty(q_list_frac.shape) its_zb = np.empty(q_list_frac.shape[0], dtype=np.int8) + count1, count2 = 0, 0 for qi,q in enumerate(q_list): - q_list_frac_fixed[qi] = np.round(CC.Methods.cart_to_cryst(rcell, _map_q_to_1st_bz(rcell, q)),6) + q_list_frac_fixed[qi] = np.round(CC.Methods.cart_to_cryst(rcell, _map_q_to_frac_bz(rcell, q)[0]),6) if np.all((np.abs(np.round(2*q_list_frac_fixed[qi], 6))%1)<1e-6): its_zb[qi] = 0 else: for alpha in range(3): - if np.abs(q_list[qi,alpha])>1e-6: - if q_list[qi,alpha] > 1e-6: + if np.abs(q_list_frac_fixed[qi,alpha])>1e-6 and np.abs(q_list_frac_fixed[qi,alpha]+0.5)>1e-6: + if q_list_frac_fixed[qi,alpha] > 1e-6: its_zb[qi] = 1 # Non zone border. Positive class. + count1 += 1 break else: its_zb[qi] = 2 # Non zone border. Negative class. + count2 += 1 break + if count1 != count2: + print("q-point classification failing...") + raise ValueError mapping = np.zeros([len(q_list), rot_cart.shape[0]], dtype=np.int32) for qi, q in enumerate(q_list): for isym in range(rot_cart.shape[0]): q_sym_cart = rot_cart[isym] @ q - q_sym_cart_1bz = _map_q_to_1st_bz(rcell, q_sym_cart) + q_sym_cart_1bz = _map_q_to_frac_bz(rcell, q_sym_cart)[0] q_sym = np.round(CC.Methods.cart_to_cryst(rcell, q_sym_cart_1bz),6) match = np.all(np.abs(q_list_frac_fixed-q_sym)<1e-3, axis=1) qii = np.where(match) @@ -336,6 +342,45 @@ def _map_q_to_1st_bz(rcell, q_cart, atol=1e-5): return tied_candidates[idx[0]] +def _map_q_to_frac_bz(rcell, q_cart, atol=1e-5): + """ + Maps a Cartesian q-vector into the [-0.5, 0.5) fractional reciprocal unit cell. + Forces +0.5 to -0.5 so opposite BZ faces map consistently. + + Parameters + ---------- + rcell : np.ndarray + Reciprocal unit cell vectors [3, 3]. + q_cart : np.ndarray + Cartesian coordinates of wave-vector [3]. + atol : float + Tolerance for floating point noise near boundary. + + Returns + ------- + q_cart_mapped : np.ndarray + Mapped Cartesian coordinates [3]. + q_frac_mapped : np.ndarray + Mapped fractional coordinates in [-0.5, 0.5) [3]. + """ + q_cart = np.array(q_cart).flatten() + rcell = np.array(rcell) + + inv_rcell = np.linalg.inv(rcell) + q_frac = q_cart @ inv_rcell + + q_frac_round = np.round(q_frac, 6) + + # Shift coordinates into the [-0.5, 0.5) range + q_frac_mapped = q_frac_round - np.floor(q_frac_round + 0.5) + + # Force values near +0.5 to strictly -0.5 + q_frac_mapped = np.where(np.abs(q_frac_mapped - 0.5) < atol, -0.5, q_frac_mapped) + + q_cart_mapped = q_frac_mapped @ rcell + + return q_cart_mapped, q_frac_mapped + def _doublet_in_list(doublet, llist, nlist): """ Return True if doublet is found in llist[:,:nlist]. diff --git a/SCHAModules/module_hess.f90 b/SCHAModules/module_hess.f90 index 77cb0df7..dadfb429 100644 --- a/SCHAModules/module_hess.f90 +++ b/SCHAModules/module_hess.f90 @@ -1,15 +1,3 @@ -! This subroutine calculates the stochastic average of an element of the -! third order force constants based on fu2 and lmat. -! The force constants element is given with three indexes that correspond -! to atom-cartesian coordinates. -! -! NOTE: there is a first addend that is not implemented at the moment -! because it should be zero if the odd correction is calculated -! at the atomic positions that minimize the gradient with respect -! to Wyckoff positions. This non-implemented addend vanishes -! if positions are fixed by symmetry. -! Version de get_v3_element_sym pero guardando datos en una matrix auxiliar - module module_hess use omp_lib @@ -138,6 +126,8 @@ subroutine get_ref3fc(nat, orbit3t, indep_3fc_elem, n_indep_3fc_elem, kernel_3fc logical, intent(in) :: verbose + double precision :: tstart, tend + double precision, dimension(nref3, 27), intent(out) :: ref_3fc integer :: nref3, dim2, nat_sc, n_mode, n_random, nr @@ -153,8 +143,9 @@ subroutine get_ref3fc(nat, orbit3t, indep_3fc_elem, n_indep_3fc_elem, kernel_3fc tmp_ref_3fc = 0 if (verbose) then + tstart = omp_get_wtime() print *, "======================= get_ref3fc() ======================= " - print *, "Computing the independent elements for the third order FCs." + print *, " Computing the independent elements for the third order FCs." end if !$omp parallel private (i, nat1, nat2, nat3, alpha, beta, gamma, v3, indep_3fc, tmp_ref_3fc) @@ -180,6 +171,9 @@ subroutine get_ref3fc(nat, orbit3t, indep_3fc_elem, n_indep_3fc_elem, kernel_3fc !$omp end do !$omp end parallel if (verbose) then + tend = omp_get_wtime() + print*, "" + print "(A, ES16.6,A)", " Elapsed time inside get_ref3fc():", tend-tstart, " seconds" print *, "======================= DONE =======================" print *, "" end if @@ -209,6 +203,8 @@ subroutine get_ref4fc(orbit4t, indep_4fc_elem, n_indep_4fc_elem, kernel_4fc, rot logical, intent(in) :: verbose double precision, dimension(nref4, 81), intent(out) :: ref_4fc + + double precision :: tstart, tend integer :: nref4, dim2, nat_sc, n_mode, n_random, nr double precision :: v4 @@ -224,6 +220,7 @@ subroutine get_ref4fc(orbit4t, indep_4fc_elem, n_indep_4fc_elem, kernel_4fc, rot tmp_ref_4fc = 0 if (verbose) then + tstart = omp_get_wtime() print *, "======================= get_ref4fc() ======================= " print *, " Computing the independent elements for the 4th order FCs" end if @@ -253,6 +250,9 @@ subroutine get_ref4fc(orbit4t, indep_4fc_elem, n_indep_4fc_elem, kernel_4fc, rot !$omp end do !$omp end parallel if (verbose) then + tend = omp_get_wtime() + print*, "" + print "(A, ES16.6,A)", " Elapsed time inside get_ref4fc():", tend-tstart, " seconds" print *, "======================= DONE =======================" print *, "" end if @@ -1081,6 +1081,8 @@ subroutine get_indep2fc( & double precision, dimension(Nqpoint, 3), intent(in) :: T_list, q_list double precision, dimension(Nqpoint, Nqpoint, n_mode, n_mode), intent(in) :: F + double precision :: tstart, tend + logical, intent(in) :: verbose double precision, dimension(nref2,9), intent(out) :: indep_fc @@ -1096,6 +1098,7 @@ subroutine get_indep2fc( & complex :: Vsa, Vsb if (verbose) then + tstart = omp_get_wtime() print*, "======================= get_indep2fc() =======================" print*, "Computing independent third order corrections to FCs \Phi^{3}_{a,b}^{\alpha,\beta}" end if @@ -1143,6 +1146,9 @@ subroutine get_indep2fc( & end do end do if (verbose) then + tend = omp_get_wtime() + print*, "" + print "(A, ES16.6,A)", " Elapsed time inside get_indep2fc():", tend-tstart, " seconds" print *, "======================= DONE =======================" print *, "" end if @@ -1413,7 +1419,7 @@ subroutine get_rotq_ws(q, pol_vecs, v4, T_list, q_list, map_tr, map_uc, isym, ro end subroutine get_rotq_ws - ! This subroutine calculates the stochastic average of an element of the +! This subroutine calculates the stochastic average of an element of the ! third order force constants based on fu2 and lmat. ! The force constants element is given with three indexes that correspond ! to atom-cartesian coordinates. @@ -1446,7 +1452,7 @@ subroutine get_v3_element_sym (na_in, nb_in, nc_in, ur, eprod, f, rho, log_err, integer :: nat_sc, n_mode, n_random, nr double precision, dimension(:), allocatable :: fun - double precision, dimension(:,:), allocatable :: Rot,v3_tp, v3_s!aux matrix where the elements have transl and perm sym applied + double precision, dimension(:,:), allocatable :: Rot, v3_tp, v3_s !aux matrix where the elements have transl and perm sym applied double precision :: v3_aux, av, av_err double precision :: aux_elem1, aux_elem2, aux_elem3 @@ -2213,6 +2219,8 @@ subroutine recognize_triplet(nat, nat_sc, tot3, nref3, nsym, mappings, map_uc, n double precision, allocatable, dimension(:,:) :: constrain_reduced + double precision :: tstart, tend + integer, dimension(tot3,3) :: all3 integer, dimension(6*nsym, 3) :: equilist integer, dimension(27) :: indep @@ -2227,7 +2235,11 @@ subroutine recognize_triplet(nat, nat_sc, tot3, nref3, nsym, mappings, map_uc, n character(len=100) :: filename - + if (verbose) then + tstart = omp_get_wtime() + print*, "=================== recognize_triplet() ===================" + print*, " Classifying atomic triplets..." + end if kernel = 0 orbit3t = 0 orbit3o = 0 @@ -2295,12 +2307,12 @@ subroutine recognize_triplet(nat, nat_sc, tot3, nref3, nsym, mappings, map_uc, n call gauss_jordan(constrain_reduced, max(nconstrain,27), 27, nconstrain, kern, indep, nindep) deallocate(constrain_reduced) - if (verbose) then - print *, "Reference triplet: (", ii, ",", jj, ",", kk, ")." - print *, "Orbit size:", equiv, ". Number of constrains: ", & -nconstrain, ". Number of independent elemements:", nindep - print *, "" - end if +! if (verbose) then +! print *, "Reference triplet: (", ii, ",", jj, ",", kk, ")." +! print *, "Orbit size:", equiv, ". Number of constrains: ", & +! nconstrain, ". Number of independent elemements:", nindep +! print *, "" +! end if do iaux = 1, 27 do jaux = 1, 27 kernel(ref3,iaux,jaux) = kern(iaux,jaux) @@ -2315,6 +2327,13 @@ subroutine recognize_triplet(nat, nat_sc, tot3, nref3, nsym, mappings, map_uc, n end do dokk end do dojj end do doii + if (verbose) then + tend = omp_get_wtime() + print "(A, I8)", " Number of independent triplets:", ref3 + print "(A, ES16.6,A)", " Elapsed time inside recognize_triplet():", tend-tstart, " seconds" + print *, "======================= DONE =======================" + print *, "" + end if end subroutine subroutine generate_rot4(rot_cart, Rot, nsym) @@ -2394,6 +2413,8 @@ subroutine recognize_quadruplet(nat, nat_sc, tot4, nref4, nsym, mappings, map_uc double precision, dimension(81,81) :: kern double precision, dimension(24*nsym*81,81) :: constrain + double precision :: tstart, tend + integer, dimension(24,4) :: permutations integer, dimension(4) :: qplet, qplet_perm, qplet_sym integer :: ii, jj, kk, ll, nall4, equiv, nconstrain, iperm, isym, iaux, & @@ -2402,7 +2423,11 @@ subroutine recognize_quadruplet(nat, nat_sc, tot4, nref4, nsym, mappings, map_uc character(len=100) :: filename - + if (verbose) then + tstart = omp_get_wtime() + print*, "================== recognize_quadruplet() ==================" + print*, " Classifying atomic quadruplets..." + end if kernel = 0 orbit4t = 0 orbit4o = 0 @@ -2477,13 +2502,13 @@ subroutine recognize_quadruplet(nat, nat_sc, tot4, nref4, nsym, mappings, map_uc call gauss_jordan(& constrain_reduced, max(nconstrain,81), 81, nconstrain, kern, indep, nindep) deallocate(constrain_reduced) - if (verbose) then - print*, & - "Reference quadruplet: (", ii, ",", jj, ",", kk, ",", ll, ")." - print*, "Orbit size:", equiv, ". Number of constrains: ", nconstrain, & - ". Number of independet elements:", nindep - print*, "" - end if + ! if (verbose) then + ! print*, & + ! "Reference quadruplet: (", ii, ",", jj, ",", kk, ",", ll, ")." + ! print*, "Orbit size:", equiv, ". Number of constrains: ", nconstrain, & + ! ". Number of independet elements:", nindep + ! print*, "" + ! end if do iaux = 1, 81 do jaux = 1, 81 kernel(ref4,iaux,jaux) = kern(iaux,jaux) @@ -2499,6 +2524,14 @@ subroutine recognize_quadruplet(nat, nat_sc, tot4, nref4, nsym, mappings, map_uc end do dokk end do dojj end do doii + if (verbose) then + tend = omp_get_wtime() + print *,"" + print "(A, I8)", " Number of independent triplets:", ref4 + print "(A, ES16.6,A)", " Elapsed time inside recognize_quadruplet():", tend-tstart, " seconds" + print *, "======================= DONE =======================" + print *, "" + end if end subroutine subroutine triplet_in_list(triplet, all3, nall3, its_in_list)