diff --git a/anharmonic/file_IO.py b/anharmonic/file_IO.py index 0f2e0e5f..d4fdf724 100644 --- a/anharmonic/file_IO.py +++ b/anharmonic/file_IO.py @@ -1395,19 +1395,6 @@ def parse_disp_fc4_yaml(filename="disp_fc4.yaml"): return new_dataset -def parse_DELTA_FORCES(disp_dataset, - filethird='FORCES_THIRD', - filesecond='FORCES_SECOND'): - forces_third = open(filethird, 'r') - forces_second = open(filesecond, 'r') - num_atom = disp_dataset['natom'] - - for disp1 in disp_dataset['first_atoms']: - second_forces = parse_force_lines(forces_second, num_atom) - for disp2 in disp1['second_atoms']: - third_forces = parse_force_lines(forces_third, num_atom) - disp2['delta_forces'] = third_forces - second_forces - def parse_DELTA_FORCES_FOURTH(disp_dataset, file4='FORCES_FOURTH', file3='FORCES_THIRD', diff --git a/anharmonic/phonon4/__init__.py b/anharmonic/phonon4/__init__.py index d8c2af83..d8302722 100644 --- a/anharmonic/phonon4/__init__.py +++ b/anharmonic/phonon4/__init__.py @@ -1,15 +1,18 @@ import numpy as np from anharmonic.phonon4.frequency_shift import FrequencyShift from phonopy.units import VaspToTHz +from phonopy.harmonic.force_constants import set_translational_invariance, set_permutation_symmetry, get_fc2 +from phonopy.structure.symmetry import Symmetry class Phono4py: def __init__(self, fc4, supercell, primitive, - mesh, + mesh=None, band_indices=None, frequency_factor_to_THz=VaspToTHz, + is_symmetry=True, is_nosym=False, symprec=1e-3, cutoff_frequency=1e-4, @@ -18,24 +21,106 @@ class Phono4py: self._fc4 = fc4 self._supercell = supercell self._primitive = primitive - self._mesh = np.intc(mesh) + + self._mesh = None + if mesh is not None: + self._mesh = np.array(mesh, dtype='intc') + if band_indices is None: self._band_indices = [ - np.arange(primitive.get_number_of_atoms() * 3)] + np.arange(primitive.get_number_of_atoms() * 3, + dtype='intc')] else: self._band_indices = band_indices self._frequency_factor_to_THz = frequency_factor_to_THz + self._is_symmetry = is_symmetry self._is_nosym = is_nosym self._symprec = symprec self._cutoff_frequency = cutoff_frequency self._log_level = log_level self._lapack_zheev_uplo = lapack_zheev_uplo - self._band_indices_flatten = np.intc( - [x for bi in self._band_indices for x in bi]) + self._band_indices_flatten = np.array( + [x for bi in self._band_indices for x in bi], + dtype='intc') self._frequency_shifts = None + + self._symmetry = Symmetry(self._supercell, + self._symprec, + self._is_symmetry) + + def get_fc2(self): + return self._fc2 + + def set_fc2(self, fc2): + self._fc2 = fc2 + + def get_fc3(self): + return self._fc3 + + def set_fc3(self, fc3): + self._fc3 = fc3 + + def get_fc4(self): + return self._fc4 + + def set_fc4(self, fc4): + self._fc4 = fc4 + + def get_symmetry(self): + return self._symmetry + + def produce_fc4(self, + forces_fc4, + displacement_dataset, + translational_symmetry_type=0, + is_permutation_symmetry=False, + is_permutation_symmetry_fc3=False, + is_permutation_symmetry_fc2=False): + file_count = 0 + for forces, disp1 in zip(forces_fc2, disp_dataset['first_atoms']): + disp1['forces'] = forces + file_count += 1 + self._fc2 = get_fc2(self._supercell, self._symmetry, disp_dataset) + if is_permutation_symmetry_fc2: + set_permutation_symmetry(self._fc2) + if translational_symmetry_type: + set_translational_invariance( + self._fc2, + translational_symmetry_type=translational_symmetry_type) + for disp1 in disp_dataset['first_atoms']: + for disp2 in disp1['second_atoms']: + disp2['forces'] = forces_fc4[file_count] + disp2['delta_forces'] = disp2['forces'] - disp1['forces'] + file_count += 1 + + self._fc3 = get_fc3( + self._supercell, + displacement_dataset, + self._fc2, + self._symmetry, + translational_symmetry_type=translational_symmetry_type, + is_permutation_symmetry=is_permutation_symmetry_fc3, + verbose=self._log_level) + + for disp1 in disp_dataset['first_atoms']: + for disp2 in disp1['second_atoms']: + for disp3 in disp1['third_atoms']: + disp3['delta_forces'] = (forces_fc4[file_count] - + disp2['forces']) + file_count += 1 + + self._fc4 = get_fc4( + self._supercell, + displacement_dataset, + self._fc3, + self._symmetry, + is_translational_symmetry=(translational_symmetry_type > 0), + is_permutation_symmetry=is_permutation_symmetry, + verbose=self._log_level) + def set_frequency_shift(self, temperatures=None): self._interaction = FrequencyShift( self._fc4, diff --git a/anharmonic/phonon4/fc4.py b/anharmonic/phonon4/fc4.py index 36f46ae6..9598207b 100644 --- a/anharmonic/phonon4/fc4.py +++ b/anharmonic/phonon4/fc4.py @@ -51,6 +51,9 @@ def get_fc4(supercell, if is_translational_symmetry: set_translational_invariance_fc4_per_index(fc4) + if is_permutation_symmetry: + set_permutation_symmetry_fc4(fc4) + return fc4 def set_translational_invariance_fc4(fc4): diff --git a/scripts/phono4py b/scripts/phono4py index b41fdc24..8b1224bf 100755 --- a/scripts/phono4py +++ b/scripts/phono4py @@ -48,9 +48,6 @@ from phonopy.units import VaspToTHz from anharmonic.phonon3.fc3 import set_permutation_symmetry_fc3, set_translational_invariance_fc3, show_drift_fc3, get_fc3 from anharmonic.file_IO import \ parse_disp_fc3_yaml, parse_disp_fc4_yaml,\ - parse_DELTA_FORCES, parse_DELTA_FORCES_FOURTH,\ - write_DELTA_FC2_SETS, parse_DELTA_FC2_SETS, \ - write_DELTA_FC2_FOURTH_SETS, parse_DELTA_FC2_FOURTH_SETS,\ write_FORCES_FOURTH, parse_FORCES_SECOND, \ read_fc4_from_hdf5, read_fc3_from_hdf5, read_fc2_from_hdf5, \ write_fc4_to_hdf5, write_fc3_to_hdf5, write_fc2_to_hdf5, \ @@ -89,7 +86,6 @@ parser.set_defaults(amplitude=None, band_indices=None, cell_poscar=None, factor=None, - fc2_fourth_sets_mode=False, forces_fourth_mode=False, grid_points=None, is_nodiag=False, @@ -103,7 +99,6 @@ parser.set_defaults(amplitude=None, log_level=None, mesh_numbers=None, primitive_axis=None, - read_fc2_fourth=False, read_fc2=False, read_fc3=False, read_fc4=False, @@ -120,10 +115,6 @@ parser.add_option("--bi", "--band_indices", dest="band_indices", parser.add_option("-c", "--cell", dest="cell_poscar", action="store", type="string", help="Read unit cell", metavar="FILE") -parser.add_option("--create_fc2_fourth", - dest="fc2_fourth_sets_mode", - action="store_true", - help="Create DELTA_FC2_FOURTH_SETS, DELTA_FC2_SETS, and fc2.hdf") parser.add_option("--cf4", "--create_f4", dest="forces_fourth_mode", action="store_true", @@ -141,10 +132,6 @@ parser.add_option("--fc2", dest="read_fc2", action="store_true", help="Read second order force constants") -parser.add_option("--fc2_fourth", - dest="read_fc2_fourth", - action="store_true", - help="Read DELTA_FC2_FOURTH_SETS, DELTA_FC2_SETS, and fc2.hdf") parser.add_option("--fc3", dest="read_fc3", action="store_true", @@ -206,13 +193,6 @@ if options.log_level is None: else: log_level = options.log_level -# Create FC2_FOURTH_SETS -if options.fc2_fourth_sets_mode: - displacements = parse_disp_fc3_yaml() - write_DELTA_FC2_FOURTH_SETS(args, displacements) - print_end() - exit(0) - # Create FC2_FOURTH_SETS if options.forces_fourth_mode: displacements = parse_disp_fc4_yaml() @@ -258,11 +238,6 @@ unitcell = read_vasp(unitcell_filename, # Supercell and Symmetry supercell = get_supercell(unitcell, settings.get_supercell_matrix()) -symmetry = Symmetry(supercell, options.symprec) - -# Log -if log_level: - print "Spacegroup: ", symmetry.get_international_table() ############################################################### # Create supercells with displacements and exit (pre-process) # @@ -299,116 +274,69 @@ else: for vec in np.dot(supercell.get_cell(), np.linalg.inv(primitive.get_cell())): print "%5.2f"*3 % tuple(vec) - # fc2 + + if options.factor is None: + factor = VaspToTHz + else: + factor = options.factor + + mesh = settings.get_mesh_numbers() + + phono4py = Phono4py(fc4, + supercell, + primitive, + mesh, + band_indices=settings.get_band_indices(), + frequency_factor_to_THz=factor, + cutoff_frequency=1e-2, + is_nosym=options.is_nosym, + symprec=options.symprec, + log_level=log_level) + + if not options.read_fc4: # Read fc4.hdf + displacements = parse_disp_fc4_yaml() + forces_fc4 = parse_FORCES_FC4(displacements) + translational_symmetry_type = is_translational_symmetry * 1 + phono4py.produce_fc4( + forces_fc4, + displacements, + translational_symmetry_type=translational_symmetry_type, + is_permutation_symmetry=is_symmetrize_fc4_r, + is_permutation_symmetry_fc3=is_symmetrize_fc3_r, + is_permutation_symmetry_fc2=is_symmetrize_fc2) + if options.read_fc2: if log_level: print "----- Read fc2 -----" sys.stdout.flush() - if os.path.exists('fc2.hdf5'): - fc2 = read_fc2_from_hdf5() - else: - print "fc2.hdf5 not found" - if log_level: - print_end() - sys.exit(0) + phono4py.set_fc2(read_fc2_from_hdf5()) else: - if log_level: - print "----- Solve fc2 -----" - sys.stdout.flush() - - disp_dataset = parse_disp_fc4_yaml() - parse_FORCES_SECOND(disp_dataset) - fc2 = get_fc2(supercell, symmetry, disp_dataset) - - if options.is_symmetrize_fc2: - set_permutation_symmetry(fc2) - - if options.is_translational_symmetry: - set_translational_invariance(fc2) - - if not options.read_fc2: if log_level: print "----- Write fc2.hdf5 -----" - write_fc2_to_hdf5(fc2) - - # fc3 - if options.read_fc3: # Read fc3.hdf5 + write_fc2_to_hdf5(phono4py.get_fc2()) + show_drift_force_constants(phono4py.get_fc2(), name='fc2') + + if options.read_fc3: if log_level: print "----- Read fc3 -----" sys.stdout.flush() - fc3 = read_fc3_from_hdf5() - if options.is_translational_symmetry: - set_translational_invariance_fc3(fc3) + phono4py.set_fc3(read_fc3_from_hdf5()) else: - if log_level: - print "----- Solve fc3 -----" - sys.stdout.flush() - - if options.read_fc2_fourth: # fc3 from DELTA_FC2_SETS - displacements = parse_disp_fc3_yaml() - parse_DELTA_FC2_SETS(displacements) - else: # fc3 from DELTA_FORCES - displacements = parse_disp_fc3_yaml() - parse_DELTA_FORCES(displacements) - - if options.is_translational_symmetry: - translational_symmetry_type = 1 - else: - translational_symmetry_type = 0 - fc3 = get_fc3( - supercell, - displacements, - fc2, - symmetry, - translational_symmetry_type=translational_symmetry_type, - is_permutation_symmetry=options.is_symmetrize_fc3_r, - verbose=log_level) - - if options.is_symmetrize_fc3_r: - if log_level: - print "----- Symmetrize fc3 by index exchange in real space -----" - set_permutation_symmetry_fc3(fc3) - - show_drift_fc3(fc3) - - if not options.read_fc3: if log_level: print "----- Write fc3.hdf5 -----" - write_fc3_to_hdf5(fc3) + write_fc3_to_hdf5(phono4py.get_fc3()) + show_drift_fc3(phono4py.get_fc3()) - # fc4 - if options.read_fc4: # Read fc4.hdf + if options.read_fc4: if log_level: print "----- Read fc4 -----" sys.stdout.flush() - fc4 = read_fc4_from_hdf5() - if options.is_translational_symmetry: - if log_level: - print "----- Impose translational invariance " - print "condition to fc4 -----" - set_translational_invariance_fc4(fc4) + phono4py.set_fc4(read_fc4_from_hdf5()) else: - if options.read_fc2_fourth: - displacements = parse_disp_fc4_yaml() - parse_DELTA_FC2_FOURTH_SETS(displacements) - else: # fc4 from FORCES_FOURTH, FORCES_THIRD and FORCES_SECOND - displacements = parse_disp_fc4_yaml() - parse_DELTA_FORCES_FOURTH(displacements) - fc4 = get_fc4( - supercell, - displacements, - fc3, - symmetry, - is_translational_symmetry=options.is_translational_symmetry, - is_permutation_symmetry=options.is_symmetrize_fc4_r, - verbose=log_level) + write_fc4_to_hdf5(phono4py.get_fc4()) + show_drift_fc4(phono4py.get_fc4()) - if options.is_symmetrize_fc4_r: - if log_level: - print "----- Symmetrize fc4 by index exchange in real space -----" - set_permutation_symmetry_fc4(fc4) - if log_level: print "(Calculating fc4 drift...)" show_drift_fc4(fc4) @@ -418,24 +346,9 @@ else: print "----- Write fc4.hdf5 -----" write_fc4_to_hdf5(fc4) - if options.factor is None: - factor = VaspToTHz - else: - factor = options.factor - - mesh = settings.get_mesh_numbers() + # Frequency shift if mesh is not None: - phono4py = Phono4py(fc4, - supercell, - primitive, - mesh, - band_indices=settings.get_band_indices(), - frequency_factor_to_THz=factor, - cutoff_frequency=1e-2, - is_nosym=options.is_nosym, - symprec=options.symprec, - log_level=log_level) phono4py.set_frequency_shift(temperatures=settings.get_temperatures()) phono4py.set_dynamical_matrix(fc2, supercell,