From a7b7994a59afded48d3f607e1d37de58608c6f12 Mon Sep 17 00:00:00 2001 From: Changho Hong Date: Tue, 1 Jun 2021 10:31:19 +0900 Subject: [PATCH 1/3] changes in convergence test to use ALGO=all if transition metal has bad convergence. In magnetic ordering, modified scheme is applied. --- src/_version.py | 2 +- src/band.py | 18 ++-- src/config.yaml | 164 ++++++++++++++++++++++++++++++++++++ src/cutoff.py | 51 +++++++++--- src/dos.py | 17 ++-- src/hse_gap.py | 20 ++--- src/kpoint.py | 82 +++++++++++++++--- src/magnetic_ordering.py | 85 ++++++++++++++----- src/mk_supercell_modi.py | 153 ++++++++++++++++++++++++++++++++++ src/module_AF.py | 174 ++++++++++++++++++++++++++++++++++++++- src/module_subr.py | 55 +++++++++++++ src/module_vasprun.py | 30 ++++++- src/relax.py | 10 ++- src/rerun_for_metal.py | 26 ++---- 14 files changed, 780 insertions(+), 107 deletions(-) create mode 100755 src/config.yaml create mode 100755 src/mk_supercell_modi.py diff --git a/src/_version.py b/src/_version.py index 50533e3..5becc17 100644 --- a/src/_version.py +++ b/src/_version.py @@ -1 +1 @@ -__version__ = "0.9.6" +__version__ = "1.0.0" diff --git a/src/band.py b/src/band.py index 66029c3..9fe0ffd 100644 --- a/src/band.py +++ b/src/band.py @@ -1,6 +1,6 @@ ########################################### -### Date: 2018-12-05 ### -### yybbyb@snu.ac.kr ### +### Date: 2020-11-05 ### +### mk01071@snu.ac.kr ### ########################################### # This is for drawing band structure and estimating band gap. import shutil, os, sys, subprocess, yaml @@ -10,7 +10,7 @@ from module_hse import * from input_conf import set_on_off from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2020-05-12' +code_data = 'Version '+__version__+'. Modified at 2019-12-17' # Set input dir = sys.argv[1] @@ -133,8 +133,11 @@ else: mag_on = check_magnet(dir+'/relax_'+pot_type,inp_yaml['magnetic_ordering']['minimum_moment']) vasprun = make_incar_for_ncl(dir_band,mag_on,kpar,npar,vasp_std,vasp_gam,vasp_ncl) - wincar(dir_band+'/INCAR',dir_band+'/INCAR',[['NSW','0'],['LCHARG','.T.'],['ALGO','Normal']],[]) + wincar(dir_band+'/INCAR',dir_band+'/INCAR',[['NSW','0'],['LCHARG','.T.']],[]) + if os.path.isfile(dir+'/relax_'+pot_type+'/CHGCAR') and os.path.getsize(dir+'/relax_'+pot_type+'/CHGCAR') > 1000: #If electronic convergence in previous step had problem + shutil.copyfile(dir+'/relax_'+pot_type+'/CHGCAR',dir_band+'/CHGCAR') + wincar(dir_band+'/INCAR',dir_band+'/INCAR',[['ICHARG','1']],[]) # VASP calculation for CHGCAR out = run_vasp(dir_band,nproc,vasprun,mpi) if out == 1: # error in vasp calculation @@ -206,14 +209,9 @@ # Drawing band structure (nkpt_for_band is used for hse band structure.) plot_band_structure(spin,[x[-1*nkpt_for_band:] for x in Band],fermi,dir_band+'/xtic.dat',dir_band+'/xlabel.dat',[inp_band['y_min'],inp_band['y_max']],dir_band) -if not os.path.isfile(inp_yaml['program']['gnuplot']): - make_amp2_log(dir_band,'If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir_band) - try: - subprocess.call([gnuplot,dir_band+'/band.in']) - except: - make_amp2_log(dir_band,'Error occured drawing figure. Please check gnuplot.') + subprocess.call([gnuplot,dir_band+'/band.in']) with open(dir_band+'/amp2.log','r') as amp2_log: with open(dir+'/amp2.log','a') as amp2_log_tot: diff --git a/src/config.yaml b/src/config.yaml new file mode 100755 index 0000000..3794516 --- /dev/null +++ b/src/config.yaml @@ -0,0 +1,164 @@ +directory: + submit: ./Submit + output: ./Output + done: ./Done + error: ./ERROR + src_path: ./src + pot_path_gga: /data/vasp4us/pot/PBE52 + pot_path_lda: /data/vasp4us/pot/LDA52 + +program: + vasp_std: ~/Initial/vasp/vasp_std + vasp_gam: /data/vasp4us/vasp_gam + vasp_ncl: /data/vasp4us/vasp_ncl + gnuplot: /gnuplot + mpi_command: mpirun + +calculation: + kp_test: 1 + encut_test: 1 + relaxation: 1 + band: 0 + density_of_states: 0 + hse_oneshot: 0 + dielectric: 0 + effective_mass: 0 + plot: 1 + magnetic_ordering: 1 + large_off: -1 + +cif2vasp: + pot_name: + GGA: + Ce: Ce_3 + LDA: +# Ge: Ge + soc_target: +# - Bi +# - Tl +# - Pb + u_value: +# La: 7.5 +# Ce: 4.5 +# All: 0.0 + magmom: 5 + max_nelm: 100 + incar: + +vasp_parallel: + npar: 2 + kpar: 1 + +convergence_test: + enconv: 0.01 + prconv: 10 + foconv: -1 + initial_kpl: 1 + max_kpl: 20 + enstart: 200 + enstep: 50 + enmax: 1000 + potential_type: GGA # Choose on potential in GGA, LDA, HSE + incar: + +relaxation: + potential_type: + - GGA +# - LDA +# - HSE + max_iteration: 10 + converged_ionic_step: 5 + pos_warning_percent: 1 + energy: -1 + pressure: 10 + force: 0.02 +# soc_on: 1 + incar: +# EDIFF: 1e-08 + +band_calculation: + kspacing_for_band: 0.005 + type_of_kpt: all # band (only for band structure plot) + y_min: 3 # set plot region, downward from VBM + y_max: 2 # set plot region, upward from CBM + potential_type: + - GGA +# - LDA +# - HSE + relax_check: 1 + incar: + KPAR: 2 + NPAR: 2 + +hybrid_oneshot: + alpha: 0.25 + potential_type: # The potential used for lattice parameter optimization + - GGA # If one variable is inserted, we use the lattice parameter and the points of VBM and CBM with that potential. +# - - HSE # If two variables are inserted, we use the lattice parameter with above potential and the points of VBM and CBM with below potential. +# - GGA + fermi_width: 0.3 + vb_dos_min: 1 + vb_dos_max: 3 + energy_width_for_extreme: 0.5 + search_space_for_extreme: 0.05 + cutoff_df_dvb: 0.3 + band_structure_correction: 1 + incar: + +density_of_states: + potential_type: + - GGA +# - LDA +# - HSE + relax_check: 1 + kp_multiplier: 2 + y_min: 3 # set plot region, downward from VBM + y_max: 2 # set plot region, upward from CBM + incar: + +dielectric: + kp_multiplier: 2 + potential_type: + - GGA +# - LDA + incar: + EDIFF: 1e-08 +# EDIFFG: 1e-04 + metal_check: 1 + relax_check: 1 + +effective_mass: + carrier_type: + - hole + - electron + pocket_search_dist: 0.02 + grid_size_for_searching: 0.02 + grid_size_for_calculation: 0.02 + temperature_for_fermi: 300 + fermi_for_cutoff: 0.99 + potential_type: + - GGA +# - LDA +# - HSE + incar: + +magnetic_ordering: + potential_type: GGA # Choose on potential in GGA, LDA, HSE (If from relax:0, it must be same to potential_type in convergence.) + minimum_moment: 0.5 + from_relax: 1 + cutoff_for_parameter: 5.0 + tolerance: 0.01 # Tolerance to define identical pair for Ising coefficient + energy_tolerance: 0.00001 # Tolerance to find lowest energy spin ordering candidates + maximum_atoms_number: -1 + maximum_pair_type: 20 + genetic_algorithm: + min_iteration: 0 + max_iteration: -1 + selection_coeff: 2 + population: + best: 15 + crossover: 30 + mutation: 30 + random: 15 + incar: +# NPAR: 2 diff --git a/src/cutoff.py b/src/cutoff.py index b6421cf..5e4763a 100644 --- a/src/cutoff.py +++ b/src/cutoff.py @@ -1,15 +1,16 @@ #################################### -# date : 2018-12-06 # +# date : 2020-11-05 # # Author : feihoon@snu.ac.kr # # Modifier : yybbyb@snu.ac.kr # +# Modifier : mk01071@snu.ac.kr # #################################### -import os, sys, subprocess, yaml +import os, sys, subprocess, yaml, shutil from module_log import * from module_vasprun import * from module_converge import * import math from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2020-05-12' +code_data = 'Version '+__version__+'. Modified at 2019-12-17' # Set input dir = sys.argv[1] @@ -105,6 +106,8 @@ incar_for_hse(now_path+'/INCAR') incar_from_yaml(now_path,inp_conv['incar']) wincar(now_path+'/INCAR',now_path+'/INCAR',[['ENCUT',str(ENCUT)],['NSW','0'],['LWAVE','F'],['LCHARG','F']],[]) + if pygrep('ALGO',dir+'/cutoff/EN'+str(ENCUT)+"/INCAR",0,0).split()[2].upper()[0] == "A" or (ENCUT-ENSTEP>=EN_recommend and os.path.isfile(dir+'/cutoff/EN'+str(ENCUT-ENSTEP)+"/INCAR") and pygrep('ALGO',dir+'/cutoff/EN'+str(ENCUT-ENSTEP)+"/INCAR",0,0).split()[2].upper()[0] == "A"): + wincar(now_path+'/INCAR',now_path+'/INCAR',[['NSW','0'],["ALGO","All"],['LCHARG','T'],['LWAVE','F']],[]) # Running vasp out = run_vasp(now_path,nproc,vasprun,mpi) if ENCUT < EN_recommend: @@ -146,9 +149,32 @@ out = electronic_step_convergence_check(now_path) if out == 2: # electronic step is not converged. (algo = normal) - make_amp2_log(dir+'/cutoff','The calculation stops but electronic step is not converged.') - print(0) - sys.exit() +# make_amp2_log(now_path,'We change ALGO normal to All.') +# wincar(now_path+'/INCAR',now_path+'/INCAR',[['ALGO','All'],['NSW','0'],['LCHARG','T'],['LWAVE','F'],['AMIX','#'],['BMIX','#'],['AMIX_MAG','#'],['BMIX_MAG','#']],[]) + make_amp2_log(now_path,'We try ALGO All with a better guess on initial charge density.') + wincar(now_path+'/INCAR',now_path+'/INCAR',[['ALGO','All'],['NSW','0'],['LCHARG','T'],['LWAVE','F'],['AMIX','#'],['BMIX','#'],['AMIX_MAG','#'],['BMIX_MAG','#']],[]) + if os.path.isfile(now_path+'/CHGCAR') and os.path.getsize(now_path+'/CHGCAR') > 1000: + wincar(now_path+'/INCAR',now_path+'/INCAR',[['ICHARG','1']],[]) + out = run_vasp(now_path,nproc,vasprun,mpi) + if out == 1: # error in vasp calculation + print(0) + sys.exit() + else: + out = electronic_step_convergence_check_CHGCAR_conv(now_path) + count = 0 + if out == 1: + count = 0; out_c = 1 + while count < 3 and out_c != 0: + out = run_vasp(now_path,nproc,vasprun,mpi) + count = count+1 + if out == 1: # error in vasp calculation + print(0) + sys.exit() + out_c = electronic_step_convergence_check_CHGCAR_conv(now_path) + if count==3 and out_c == 1: + make_amp2_log(target,'It is not converged.') + print(0) + sys.exit() write_conv_result(now_path,enlog) # electronic step is converged. @@ -158,6 +184,12 @@ os.chdir('../') make_amp2_log(dir+'/cutoff','cut off energy test is done') +#cp converged CHGCAR to INPUT0 +if os.path.isfile(dir+'/cutoff/EN'+str(ENCUT-3*ENSTEP)+'/CHGCAR') and os.path.getsize(dir+'/cutoff/EN'+str(ENCUT-3*ENSTEP)+'/CHGCAR') > 1000: + shutil.copyfile(dir+'/cutoff/EN'+str(ENCUT-3*ENSTEP)+'/CHGCAR',dir+'/INPUT0/CHGCAR_conv') +if pygrep('ALGO',dir+'/cutoff/EN'+str(ENCUT-3*ENSTEP)+"/INCAR",0,0).split()[2].upper()[0] == "A": + wincar(dir+'/INPUT0/INCAR',dir+'/INPUT0/INCAR',["ALGO","All"],[]) + with open(enlog,'a') as result: result.write("\nConvergence criterion: E/atom < "+str(ENCONV)+" eV\n") @@ -169,14 +201,9 @@ wincar(dir+'/INPUT0/INCAR',dir+'/INPUT0/INCAR',[['ENCUT',str(ENCUT-3*ENSTEP)]],[]) make_conv_dat(dir+'/cutoff','cutoff') -if not os.path.isfile(inp_yaml['program']['gnuplot']): - make_amp2_log(dir+'/cutoff','If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir+'/cutoff') - try: - subprocess.call([gnuplot,dir+'/cutoff/conv_plot.in']) - except: - make_amp2_log(dir+'/cutoff','Error occured drawing figure. Please check gnuplot.') + subprocess.call([gnuplot,dir+'/cutoff/conv_plot.in']) with open(dir+'/cutoff/amp2.log','r') as amp2_log: with open(dir+'/amp2.log','a') as amp2_log_tot: diff --git a/src/dos.py b/src/dos.py index 03ec766..8d70f31 100644 --- a/src/dos.py +++ b/src/dos.py @@ -1,6 +1,6 @@ #################################### -# Modifier : yybbyb@snu.ac.kr # -# data : 2018-12-05 # +# Date : 2020-11-05 # +# mk01071@snu.ac.kr # #################################### # This is for drawing density of states. import shutil, os, sys, subprocess, yaml @@ -9,7 +9,7 @@ from module_dos import * from input_conf import set_on_off from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2020-05-12' +code_data = 'Version '+__version__+'. Modified at 2019-12-17' # Set input dir = sys.argv[1] @@ -176,7 +176,7 @@ sys.exit() make_amp2_log(dir_dos,'CHGCAR file is generated successfully.') - wincar(dir_dos+'/INCAR',dir_dos+'/INCAR',[['NEDOS','3001'],['ISMEAR','-5'],['SIGMA','0.1']],[]) + wincar(dir_dos+'/INCAR',dir_dos+'/INCAR',[['NEDOS','3001'],['ISMEAR','-5'],['SIGMA','0.1'],["ALGO","Normal"]],[]) incar_from_yaml(dir_dos,inp_dos['incar']) @@ -189,7 +189,7 @@ while 1: make_multiple_kpts(dir+'/kptest/kpoint.log',dir_dos+'/KPOINTS',dir_dos+'/POSCAR',inp_dos['kp_multiplier'],sym,gam_option) - wincar(dir_dos+'/INCAR',dir_dos+'/INCAR',[['NSW','0'],['ISTART','1'],['ICHARG','11'],['LCHARG','.F.']],[]) + wincar(dir_dos+'/INCAR',dir_dos+'/INCAR',[['NSW','0'],['ICHARG','11'],['LCHARG','.F.']],[]) if no_rlx == 1: mag_on = 2 else: @@ -247,14 +247,9 @@ make_dos_in(dir_dos,atom_name,spin,len(par_dos[0][0]),[inp_dos['y_min'],inp_dos['y_max']+gap]) -if not os.path.isfile(inp_yaml['program']['gnuplot']): - make_amp2_log(dir_dos,'If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir_dos+'/Pdos_dat') - try: - subprocess.call([gnuplot,dir_dos+'/Pdos_dat/dos.in']) - except: - make_amp2_log(dir_dos,'Error occured drawing figure. Please check gnuplot.') + subprocess.call([gnuplot,dir_dos+'/Pdos_dat/dos.in']) make_amp2_log(dir_dos,'DOS calculation is done.') diff --git a/src/hse_gap.py b/src/hse_gap.py index cbd9c97..71a209c 100644 --- a/src/hse_gap.py +++ b/src/hse_gap.py @@ -12,7 +12,7 @@ from input_conf import set_on_off from module_relax import * from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2020-05-12' +code_data = 'Version '+__version__+'. Modified at 2019-12-17' # Set input dir = sys.argv[1] @@ -171,9 +171,9 @@ incar_for_hse(dir_hse+'/INCAR') if PBE0_on == 1: inp_hse['alpha'] = calc_alpha_auto(diel_path+'/dielectric.log') - wincar(dir_hse+'/INCAR',dir_hse+'/INCAR',[['NSW','0'],['HFSCREEN','0.0'],['ALGO','All'],['AEXX',str(inp_hse['alpha'])]],[]) + wincar(dir_hse+'/INCAR',dir_hse+'/INCAR',[['NSW','0'],['HFSCREEN','0.0'],['ALGO','ALL'],['AEXX',str(inp_hse['alpha'])]],[]) else: - wincar(dir_hse+'/INCAR',dir_hse+'/INCAR',[['NSW','0'],['HFSCREEN','0.2'],['ALGO','All'],['AEXX',str(inp_hse['alpha'])]],[]) + wincar(dir_hse+'/INCAR',dir_hse+'/INCAR',[['NSW','0'],['HFSCREEN','0.2'],['ALGO','ALL'],['AEXX',str(inp_hse['alpha'])]],[]) incar_from_yaml(dir_hse,inp_hse['incar']) mag_on = check_magnet(dir+'/relax_'+pot_cell,inp_yaml['magnetic_ordering']['minimum_moment']) vasprun = make_incar_for_ncl(dir_hse,mag_on,kpar,npar,vasp_std,vasp_gam,vasp_ncl) @@ -219,14 +219,9 @@ for k in range(len(KPT)): Band[n][k][i] = Band[n][k][i] + E_shift plot_band_corrected_structure(spin,Band,eVBM,dir_band+'/xtic.dat',dir_band+'/xlabel.dat',[inp_band['y_min'],inp_band['y_max']+E_shift],dir_band) - if not os.path.isfile(inp_yaml['program']['gnuplot']): - make_amp2_log(dir_hse,'If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir_band) - try: - subprocess.call([gnuplot,dir_band+'/band_corrected.in']) - except: - make_amp2_log(dir_hse,'Error occured drawing figure. Please check gnuplot.') + subprocess.call([gnuplot,dir_band+'/band_corrected.in']) # Metal in GGA. Band reordering is required. else: @@ -257,14 +252,9 @@ fermi = get_fermi_level(Band_reorder,nelect,ncl) if calc_gap(fermi,spin,ncl,KPT[:num_kpt_for_image],Band_reorder,nelect) > 0.01: plot_band_corrected_structure(spin,Band_reorder,eVBM,dir_band+'/xtic.dat',dir_band+'/xlabel.dat',[inp_band['y_min'],inp_band['y_max']+float(gap)],dir_band) - if not os.path.isfile(inp_yaml['program']['gnuplot']): - make_amp2_log(dir_hse,'If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir_band) - try: - subprocess.call([gnuplot,dir_band+'/band_corrected.in']) - except: - make_amp2_log(dir_hse,'Error occured drawing figure. Please check gnuplot.') + subprocess.call([gnuplot,dir_band+'/band_corrected.in']) else: make_amp2_log(dir_hse,'Warning. We cannot correct the band structure.') else: diff --git a/src/kpoint.py b/src/kpoint.py index a041e32..88bc8e9 100644 --- a/src/kpoint.py +++ b/src/kpoint.py @@ -1,14 +1,15 @@ #################################### -# date : 2018-12-06 # +# Date: 2020-11-05 # # Author : feihoon@snu.ac.kr # # Modifier : yybbyb@snu.ac.kr # +# Modifier : mk01071@snu.ac.kr # #################################### import os, sys, subprocess, yaml from module_log import * from module_vasprun import * from module_converge import * from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2020-05-12' +code_data = 'Version '+__version__+'. Modified at 2019-12-17' # Set input dir = sys.argv[1] @@ -98,6 +99,10 @@ incar_for_hse(now_path+'/INCAR') incar_from_yaml(now_path,inp_conv['incar']) wincar(now_path+'/INCAR',now_path+'/INCAR',[['NSW','0'],['LCHARG','F'],['LWAVE','F']],[]) + #if algo is all in the last step, current step uses ALGO all + if os.path.isfile(dir+'/kptest/KP'+str(KPL-1)+"/INCAR") and pygrep('ALGO',dir+'/kptest/KP'+str(KPL-1)+"/INCAR",0,0).split()[2].upper()[0] == "A": + wincar(now_path+'/INCAR',now_path+'/INCAR',[['NSW','0'],["ALGO","All"],['LCHARG','T'],['LWAVE','F']],[]) + out = run_vasp(now_path,nproc,vasprun,mpi) if KPL == 1 or KPL == 2: if out == 1: # error in vasp calculation @@ -118,8 +123,38 @@ out = electronic_step_convergence_check(now_path) if out == 2: # electronic step is not converged. (algo = normal) - make_amp2_log(dir+'/kptest','The calculation stops but electronic step is not converged, but pass because of too small kpoints.') - loopnum = 0 +# make_amp2_log(dir+'/kptest','The calculation stops but electronic step is not converged, but pass because of too small kpoints.') + make_amp2_log(now_path,'We try ALGO All with a better guess on initial charge density.') + wincar(now_path+'/INCAR',now_path+'/INCAR',[['ALGO','All'],['NSW','0'],['LCHARG','T'],['LWAVE','F'],['AMIX','#'],['BMIX','#'],['AMIX_MAG','#'],['BMIX_MAG','#']],[]) + if os.path.isfile(now_path+'/CHGCAR') and os.path.getsize(now_path+'/CHGCAR') > 1000: + wincar(now_path+'/INCAR',now_path+'/INCAR',[['ICHARG','1']],[]) + out = run_vasp(now_path,nproc,vasprun,mpi) + if out == 1: # error in vasp calculation + make_amp2_log(dir+'/kptest','VASP error occurs, but pass because of too small kpoints.') + out = 3 + loopnum = 0 + else: + out = electronic_step_convergence_check_CHGCAR_conv(now_path) + count = 0 + if out == 1: + count = 0; out_c = 1 + while count < 3 and out_c != 0: + make_amp2_log(now_path,'The last CHGCAR is used as an initial guess.') + out = run_vasp(now_path,nproc,vasprun,mpi) + count = count+1 + if out == 1: # error in vasp calculation + make_amp2_log(dir+'/kptest','VASP error occurs, but pass because of too small kpoints.') + out = 3 + loopnum = 0 + break + out_c = electronic_step_convergence_check_CHGCAR_conv(now_path) + if count==3 and out_c == 1: + make_amp2_log(dir+'/kptest','The calculation stops but electronic step is not converged, but pass because of too small kpoints.') + out = 3 + loopnum = 0 + continue + write_conv_result(now_path,kplog) + elif out < 2: write_conv_result(now_path,kplog) @@ -137,9 +172,33 @@ out = electronic_step_convergence_check(now_path) if out == 2: # electronic step is not converged. (algo = normal) - make_amp2_log(dir+'/kptest','The calculation stops but electronic step is not converged.') - print(0) - sys.exit() + make_amp2_log(now_path,'We try ALGO All with a better guess on initial charge density.') + wincar(now_path+'/INCAR',now_path+'/INCAR',[['ALGO','All'],['NSW','0'],['LCHARG','T'],['LWAVE','F'],['AMIX','#'],['BMIX','#'],['AMIX_MAG','#'],['BMIX_MAG','#']],[]) + if os.path.isfile(now_path+'/CHGCAR') and os.path.getsize(now_path+'/CHGCAR') > 1000: + wincar(now_path+'/INCAR',now_path+'/INCAR',[['ICHARG','1']],[]) + out = run_vasp(now_path,nproc,vasprun,mpi) + if out == 1: # error in vasp calculation + print(0) + sys.exit() + else: + out = electronic_step_convergence_check_CHGCAR_conv(now_path) + count = 0 + if out == 1: + count = 0; out_c = 1 + while count < 3 and out_c != 0: + make_amp2_log(now_path,'The last CHGCAR is used as an initial guess.') + out = run_vasp(now_path,nproc,vasprun,mpi) + count = count+1 + if out == 1: # error in vasp calculation + print(0) + sys.exit() + out_c = electronic_step_convergence_check_CHGCAR_conv(now_path) + if count==3 and out_c == 1: + make_amp2_log(target,'It is not converged.') + print(0) + sys.exit() + os.remove(now_path+'/CHGCAR') + os.remove(now_path+'/CHG') # electronic step is converged. write_conv_result(now_path,kplog) @@ -159,16 +218,13 @@ result.write("Converged KPL: "+str(KPL-3)+'\n') subprocess.call(['cp',dir+'/kptest/KP'+str(KPL-3)+'/KPOINTS',dir+'/kptest/KPOINTS_converged']) subprocess.call(['cp',dir+'/kptest/KPOINTS_converged',dir+'/INPUT0/KPOINTS']) +if pygrep('ALGO',dir+'/kptest/KP'+str(KPL-3)+"/INCAR",0,0).split()[2].upper()[0] == "A": + wincar(dir+'/INPUT0/INCAR',dir+'/INPUT0/INCAR',[["ALGO","All"]],[]) make_conv_dat(dir+'/kptest','kpoint') -if not os.path.isfile(inp_yaml['program']['gnuplot']): - make_amp2_log(dir+'/kptest','If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir+'/kptest') - try: - subprocess.call([gnuplot,dir+'/kptest/conv_plot.in']) - except: - make_amp2_log(dir+'/kptest','Error occured drawing figure. Please check gnuplot.') + subprocess.call([gnuplot,dir+'/kptest/conv_plot.in']) with open(dir+'/kptest/amp2.log','r') as amp2_log: with open(dir+'/amp2.log','a') as amp2_log_tot: diff --git a/src/magnetic_ordering.py b/src/magnetic_ordering.py index b0ae48d..b0dc5c8 100644 --- a/src/magnetic_ordering.py +++ b/src/magnetic_ordering.py @@ -1,6 +1,6 @@ #################################### -# date : 2019-04-08 # -# Author : yybbyb@snu.ac.kr # +# date : 2020-11-05 # +# Author : kstgood333@snu.ac.kr # #################################### # This is for identifying the most stable magnetic spin ordering. @@ -36,6 +36,12 @@ inp_rlx = inp_yaml['relaxation'] cutoff_length = inp_af['cutoff_for_parameter'] pot_type = inp_af['potential_type'] +#mag_ver = inp_af['ising_model'] +#if mag_ver == 'Modified': +# MAG_V = '1' +#else: +# MAG_V = '0' +MAG_V = 1 if pot_type == 'LDA': POT = 'LDA' else: @@ -119,17 +125,38 @@ # Make supercell -if not os.path.isfile(target+'/POSCAR_param'): - min_cell_length = cutoff_length*2 - while 1: - try: - out = int(subprocess.check_output([pypath,src_path+'/mk_supercell.py',inp_pos,dir+'/INPUT0/KPOINTS',str(min_cell_length),target+'/POSCAR_param'],universal_newlines=True).split()[0]) - except: - out = 1 - if out == 0: - break - else: - min_cell_length = min_cell_length + 2 +if MAG_V == '1': + if not os.path.isfile(target+'/POSCAR_param'): + with open(target+'/tmp_mag_list.dat','w') as tml: + tml.write(' '.join(mag_atom_list)) + len_iter = 0.0 + while 1: + try: + out = int(subprocess.check_output([pypath,src_path+'/mk_supercell_modi.py',inp_pos,target+'/POSCAR_param',str(len_iter)],universal_newlines=True).split()[0]) + except: + out = 1 + if out == 0 and len_iter < 10: + break + if out == 1 and len_iter < 10: + len_iter += 0.5 + if len_iter >= 10: + make_amp2_log(target,'fail to find minimum supercell.') + print(1) + sys.exit() + os.remove(target+'/tmp_mag_list.dat') + mk_kpoints_modi(dir) +else: + if not os.path.isfile(target+'/POSCAR_param'): + min_cell_length = cutoff_length*2 + while 1: + try: + out = int(subprocess.check_output([pypath,src_path+'/mk_supercell.py',inp_pos,dir+'/INPUT0/KPOINTS',str(min_cell_length),target+'/POSCAR_param'],universal_newlines=True).split()[0]) + except: + out = 1 + if out == 0: + break + else: + min_cell_length = min_cell_length + 2 # Check whether the genetic alrogirhm starts or not [axis_param,pos_param] = read_poscar(target+'/POSCAR_param') @@ -162,8 +189,11 @@ os.chdir(targ_dir) copy_input_no_kp(dir+'/INPUT0',targ_dir,POT) subprocess.call(['cp',target+'/POSCAR_param',targ_dir+'/POSCAR']) - subprocess.call(['cp',src_path+'/KPOINTS_gamma',targ_dir+'/KPOINTS']) - wincar(targ_dir+'/INCAR',targ_dir+'/INCAR',[['LCHARG','F'],['LWAVE','F'],['ALGO','Normal'],['NELM',max_nelm]],[]) + if MAG_V == '1': + subprocess.call(['cp',dir+'/INPUT0/KPOINTS_modi',targ_dir+'/KPOINTS']) + else: + subprocess.call(['cp',src_path+'/KPOINTS_gamma',targ_dir+'/KPOINTS']) + wincar(targ_dir+'/INCAR',targ_dir+'/INCAR',[['LCHARG','F'],['LWAVE','F'],['ALGO','All'],['NELM',max_nelm]],[]) if pot_type == 'HSE': incar_for_hse(targ_dir+'/INCAR') wincar(targ_dir+'/INCAR',targ_dir+'/INCAR',[['ALGO','All']],[]) @@ -177,7 +207,7 @@ vasprun = vasp_std # VASP calculation for CHGCAR - wincar(targ_dir+'/INCAR',targ_dir+'/INCAR',[['NSW','0'],['ISYM','0'],['NELM','50'],['LCHARG','T']],[]) + wincar(targ_dir+'/INCAR',targ_dir+'/INCAR',[['NSW','0'],['ISYM','0'],['NELM','100'],['LCHARG','T']],[]) out = run_vasp(targ_dir,nproc,vasprun,mpi) if out == 1: # error in vasp calculation @@ -204,6 +234,9 @@ os.remove(targ_dir+'/CHGCAR') os.remove(targ_dir+'/CHG') + final_mag = mk_mag_list_modi(targ_dir) + with open(targ_dir+'/mag','w') as magf: + magf.write(final_mag+'\n') energy = pygrep('free ','OUTCAR',0,0).splitlines()[-1].split()[4] with open(targ_dir+'/energy','w') as enef: enef.write('Spin_'+str(i)+'\t'+energy+'\n') @@ -212,14 +245,24 @@ with open(target+'/energy','w') as enef: with open(targ_dir+'/energy','r') as ener: enef.write(ener.read()) + with open(target+'/mag','w') as magf: + with open(targ_dir+'/mag','r') as magr: + magf.write(magr.read()) else: with open(target+'/energy','a') as enef: with open(targ_dir+'/energy','r') as ener: enef.write(ener.read()) - -JJ = calc_ising_param(sole_list,pair_list,target+'/energy') -write_ising_param(JJ,target) -write_inp_for_GA(mag_atom_list,inp_af,target) + with open(target+'/mag','a') as magf: + with open(targ_dir+'/mag','r') as magr: + magf.write(magr.read()) +if MAG_V == '1': + JJ = calc_ising_param_modi(mag_atom_list,target) + write_ising_param_modi(JJ,target) + write_inp_for_GA(mag_atom_list,inp_af,target) +else: + JJ = calc_ising_param(sole_list,pair_list,target+'/energy') + write_ising_param(JJ,target) + write_inp_for_GA(mag_atom_list,inp_af,target) ### genetic algorithm in supercell supercell_list = set_supercell_list(inp_pos,mag_atom_list) @@ -258,7 +301,7 @@ make_amp2_log(target,'Smallest cell is in GA_'+''.join([str(x) for x in optimized_supercell])+'.') ground_path = target+'/GA_'+''.join([str(x) for x in optimized_supercell]) os.chdir(ground_path) -POSCARs = glob.glob(ground_path+'/POSCAR_spin*') +POSCARs = sorted(glob.glob(ground_path+'/POSCAR_spin*')) pos_atom_num = [] for i in range(len(POSCARs)): if not os.path.isdir(ground_path+'/Stable'+str(i)): diff --git a/src/mk_supercell_modi.py b/src/mk_supercell_modi.py new file mode 100755 index 0000000..c7ad299 --- /dev/null +++ b/src/mk_supercell_modi.py @@ -0,0 +1,153 @@ +# This is a code to build supercell for the Ising coefficient +import os, math, sys, yaml +import numpy as np +import module_subr +from itertools import permutations +from itertools import combinations +import operator +from time import strftime +t_start = strftime("%y%m%d-%H%M%S") +try: + pos_file = str(sys.argv[1]) +except: + sys.exit() + +len_min = 5.0 +out_pos_file = sys.argv[2] +len_iter = float(sys.argv[3]) +latt_diff = 1.0 + +len_crit = len_min+len_iter+latt_diff +min_natom = 16 # number of atoms in primitive cell +max_angle = 100 +min_angle = 60 +# Using relaxed position file +poslines = module_subr.read_file(pos_file) + +# Read lattice vector & atom position +lat_prim=module_subr.read_lat(poslines) +atnames= poslines[5].split() +totnum= sum([int(x) for x in poslines[6].split()]) +at_prim = module_subr.read_at_pos(poslines,totnum) +vol_prim = np.inner(lat_prim[0],np.cross(lat_prim[1],lat_prim[2])) + +# Make combination of rotation vector ex) [-2,1,2] +m = 1 +len_check = 1 +tmp = [] +bf_list = [] +# only min_length check +while len_check: + rot = module_subr.mk_combination(m) + if m!=1: + rot_small = module_subr.mk_combination(m-1) + rot_new = [x for x in rot if x not in rot_small] + bf_list = new_list + rot = rot_new + + # Remove duplicate or all element < 0 + rot_sort = rot + # Check length criterion + [new_list,min_sw] = module_subr.check_length(rot_sort, lat_prim, len_min, len_crit) + if (len(bf_list) !=0) and (len(new_list) == 0): #until do not append new list, it takes loop infinitely + len_check = 0 + else: + if min_sw == 0: + print(1) + sys.exit() + m = m+1 + + tmp = tmp+new_list + + +new_list = tmp +# make rotation matrix using basis +# format : [[rot_mat],[lattice_mat]] +basis_mat = module_subr.mk_rot_matrix(new_list) + +# remove det = 0 +basis_mat = module_subr.rm_det_0(basis_mat) +# remove too sharp shape structures +basis_mat = module_subr.check_shape(basis_mat, min_angle, max_angle) + +# Check the number of atoms: +basis_mat = module_subr.check_natom(basis_mat, lat_prim, totnum, min_natom) + +# Check defect distance: +basis_mat = module_subr.check_defect_length(basis_mat, len_min) + +# Check pair: +#cri_natom = 100 +#best_set = [] +#lattice diff check and tot number check +#for i in range(len(basis_mat)): +# lat_mat = basis_mat[i][1] +# length = module_subr.cal_lat_len(lat_mat) +# tot_atom = module_subr.cal_volume(lat_mat)/vol_prim*totnum + +# if (max(length)-min(length))/max(length) < latt_diff: +# if tot_atom < cri_natom : +# best_set.append([basis_mat[i][0],basis_mat[i][1]]) + + +#find small natoms +#basis_mat = best_set +target = sys.argv[2].split('POSCAR_param')[0] +with open (target+'tmp_mag_list.dat','r') as f1: + mag_atom_list = f1.readlines()[0].split() + +check_set = module_subr.check_pair(basis_mat,at_prim,mag_atom_list) + +best_natom = 10e5 +best = 'Null' +best_set = [] + + +for i in range(len(check_set)): + lat_mat = check_set[i][1] + length = module_subr.cal_lat_len(lat_mat) + tot_atom = module_subr.cal_volume(lat_mat)/vol_prim*totnum + + if tot_atom < best_natom: + best_natom = tot_atom + best_set = [[check_set[i][0],check_set[i][1]]] + elif tot_atom == best_natom: + best_set.append([check_set[i][0],check_set[i][1]]) +#best_set = check_set + + +#best_set = module_subr.check_pair(best_set,at_prim,mag_atom_list) +best_latt_sum = -10e5 + +for i in range(len(best_set)): + latt_sum = 0 + for j in range(3): + for k in range(3): + latt_sum = latt_sum + best_set[i][1][j][k] + + if latt_sum > best_latt_sum: + best_latt_sum = latt_sum + best = best_set[i] +det = int(round(best_natom))//int(totnum) +# Set atomic position +if best == 'Null': + print(1) + sys.exit() +else: + at = module_subr.set_atom(at_prim, best[0]) +lat_mat = best[1] +# Create POSCAR files (GGA/HSE) +with open(out_pos_file,'w') as f1: + f1.write('Clean supercell of '+str(atnames)+'\n') + f1.write(' 1.00000000000000\n') + f1.write("{:21.13f}{:19.13f}{:19.13f}\n".format(lat_mat[0][0],lat_mat[0][1],lat_mat[0][2])) + f1.write("{:21.13f}{:19.13f}{:19.13f}\n".format(lat_mat[1][0],lat_mat[1][1],lat_mat[1][2])) + f1.write("{:21.13f}{:19.13f}{:19.13f}\n".format(lat_mat[2][0],lat_mat[2][1],lat_mat[2][2])) + f1.write(poslines[5]) + f1.write(' '+' '.join([str(int(x)*det) for x in poslines[6].split()])+'\n') + f1.write('Selective dynamics\nDirect\n') + for t in range(len(at)): + f1.write("{:19.15f}{:19.15f}{:19.15f} T T T\t! {}\n".format(at[t][0],at[t][1],at[t][2],at[t][3])) + +t_end = strftime("%y%m%d-%H%M%S") +print(0) diff --git a/src/module_AF.py b/src/module_AF.py index d291cdb..c163725 100644 --- a/src/module_AF.py +++ b/src/module_AF.py @@ -3,11 +3,12 @@ # Author : yybbyb@snu.ac.kr # #################################### # This is a package of modules for identifying the most stable magnetic spin ordering. -import os,sys,yaml,subprocess +import os,sys,yaml,subprocess, re import numpy as np from module_vector import * from module_vasprun import * from module_amp2_input import * +from module_GA import * def check_spin(ref_dir,inp_pos,min_mom): [axis,atom_pos] = read_poscar(inp_pos) @@ -103,6 +104,86 @@ def reduce_mag_list(mag_list): spin.append(i) num_spin.append(1) return [spin,num_spin] +## for modified Ising model ## +def mk_p_list(target): + with open(target+'/pair_list_table.dat','r') as f1: + plt_lines = f1.readlines() + p_list = [] + for i in range(len(plt_lines)): + tmp_list = [0 for j in range(3)] + tmp = plt_lines[i].split() + for j in range(1,3): + tmp_list[j-1] = float(re.findall('\d+',tmp[j])[0]) + tmp_list[2] = round(float(tmp[0]),2) + p_list.append(tmp_list[0:3]) + return p_list + +def mk_modi_ising_list(mag_lines,p_list,target,mag_atom_list,axis,atom_pos): + pattern = re.compile(r'((-)?\d{1,3}(,\d{3})*(\.\d+)?)') + m_atoms = mag_lines.split() + m_atoms = list(map(float,[x.split('*')[0] for x in mag_lines.split()])) + s_atoms = list(map(float,[x.split('*')[1] for x in mag_lines.split()])) + atom_pos_cart = [] + for j in range(len(m_atoms)): + for k in range(int(sum(m_atoms[0:j])),int(sum(m_atoms[0:j+1]))): + if s_atoms[j] < 0: + atom_pos_cart.append(dir_to_cart(atom_pos[k][0:3],axis)+atom_pos[k][3:]) + atom_pos_cart[k][-1] = atom_pos_cart[k][-1]+'_down' + else: + atom_pos_cart.append(dir_to_cart(atom_pos[k][0:3],axis)+atom_pos[k][3:]) + atom_pos_cart[k][-1] = atom_pos_cart[k][-1]+'_up' + mag_true = [] + mag_sum = 0 + for j in range(len(atom_pos_cart)): + if atom_pos_cart[j][4] in mag_atom_list: + mag_true.append(1) + mag_sum = mag_sum +1 + else: + mag_true.append(0) + info_pair = pairCategorize(atom_pos_cart,axis,mag_true,5.0) + count = {} + for line in info_pair: + if line[5] in list(count.keys()): + count[line[5]] = count[line[5]] +1 + else: + count[line[5]] = 1 + keys = list(count.keys()) + tot_list=[] + for j in range(len(keys)): #info_pair to list, format ex) [[1,2,3.66(pair length),4(count),Fe1_upFe2_down3.66(info_pair key)]] + tmp_list = [] + idx_check = pattern.findall(keys[j]) + if idx_check[0][0] >= idx_check[1][0]: + tmp_list.append(float(idx_check[1][0])) + tmp_list.append(float(idx_check[0][0])) + tmp_list.append(float(idx_check[2][0])) + else: + tmp_list.append(float(idx_check[0][0])) + tmp_list.append(float(idx_check[1][0])) + tmp_list.append(float(idx_check[2][0])) + tmp_list.append(count[keys[j]]) + tmp_list.append(keys[j]) + tot_list.append(tmp_list) + + m_list = [[] for j in range(len(p_list))] + for j in range(len(tot_list)): + for k in range(len(p_list)): + if tot_list[j][0:3] == p_list[k]: + m_list[k].append(tot_list[j]) + return m_list + +def cal_pinv(e_list,N_list,target): + e_array = array(e_list) + N_array = array(N_list) + JJ = dot(linalg.pinv(N_array),e_array) + return JJ + +def mk_e_list(target): + with open(target+'/energy','r') as f1: + e_lines = f1.readlines() + e_list = [] + for i in range(len(e_lines)): + e_list.append(float(e_lines[i].split()[-1])) + return e_list # This function is for calculating ising parameter from magnetic pairs def calc_ising_param(sole_list,pair_list,ene_file): @@ -124,6 +205,48 @@ def write_ising_param(JJ,target): for line in JJ: pcw.write(str(line[0])+'\t'+str(line[3])+'\t'+line[1]+'\t'+line[2]+'\n') +# This function is for calculating modified ising parameter from magnetic pairs +def calc_ising_param_modi(mag_atom_list,target): + up_p = re.compile(r'up') + down_p = re.compile(r'down') + spin_atom_list = [] + for i in range(len(mag_atom_list)): + spin_atom_list.append(mag_atom_list[i]+'_up') + spin_atom_list.append(mag_atom_list[i]+'_down') + p_list = mk_p_list(target) + [axis, atom_pos] = read_poscar(target+'/POSCAR_param') + with open(target+'/mag','r') as f1: + mag_lines = f1.readlines() + N_list = [] + for i in range(len(mag_lines)): + # m_list format = [[1,1,3.82,4,Fe1_upFe1_up3.82][1,1,3.02,6,Fe1_upFe1_down3.02]]..[mag_atom_a,mag_atom_b,pair_length,count,spin_info] + m_list = mk_modi_ising_list(mag_lines[i],p_list,target,spin_atom_list,axis,atom_pos) + N_tmp = [1] + for j in range(len(p_list)): + ferro = 0 + af = 0 + for k in range(len(m_list[j])): + a = up_p.findall(m_list[j][k][-1]) + b = down_p.findall(m_list[j][k][-1]) + if len(a) == 2 or len(b) == 2: + ferro += m_list[j][k][-2] + if len(a) == 1 and len(b) == 1: + af += m_list[j][k][-2] + tmp = (ferro-af)/2 + N_tmp.append(tmp) + N_list.append(N_tmp) + e_list = mk_e_list(target) + JJ = cal_pinv(e_list,N_list,target) + return JJ + +def write_ising_param_modi(JJ,target): + with open(target+'/pair_list_table.dat','r') as f1: + pt_lines = f1.readlines() + with open(target+'/Pair_coeff.dat','w') as f1: + for i in range(len(pt_lines)): + p_lines = pt_lines[i].split() + f1.write('{}\t {} \t {} \t {} \n'.format(p_lines[0],JJ[i+1],p_lines[1],p_lines[2])) + # This function is for writing inputs for genetic algorithm def write_inp_for_GA(mag_atom_list,inp_af,target): inp_yaml = {} @@ -224,3 +347,52 @@ def set_mag_info(atom_pos,target,mag_val): [spin,num_spin] = reduce_mag_list(magmom) with open(target+'/mag_info','w') as inp: inp.write(' '.join([str(num_spin[x])+'*'+str(spin[x]) for x in range(len(spin))])) + +def mk_kpoints_modi(c_dir): + axis_conv = poscar_to_axis(c_dir+'/relax_GGA/POSCAR') + recipro_latt_conv = reciprocal_lattice(axis_conv) + l_conv = [] + for i in range(3): + l_conv.append((recipro_latt_conv[i][0]**2.+recipro_latt_conv[i][1]**2.+recipro_latt_conv[i][2]**2.)**0.5) + axis_modi = poscar_to_axis(c_dir+'/magnetic_ordering/POSCAR_param') + recipro_latt_modi = reciprocal_lattice(axis_modi) + l_modi = [] + for i in range(3): + l_modi.append((recipro_latt_modi[i][0]**2.+recipro_latt_modi[i][1]**2.+recipro_latt_modi[i][2]**2.)**0.5) + l_ratio = [l_modi[x]/l_conv[x] for x in range(3) ] + with open(c_dir+'/INPUT0/KPOINTS','r') as f1: + kp_lines = f1.readlines() + kp_conv = kp_lines[3].split() + kp_conv = list(map(int,kp_conv)) + with open(c_dir+'/INPUT0/KPOINTS_modi','w') as kp_modi: + for j in range(3): + kp_modi.write(kp_lines[j]) + KP = [] + for j in range(3): + KP.append(str(int(round(kp_conv[j]*l_ratio[j])))) + if KP[-1] == '0': + KP[-1] = '1' + kp_modi.write(" "+KP[0]+" "+ KP[1]+" "+KP[2]+"\n 0 0 0\n") + +def mk_mag_list_modi(targ_dir): + datfile = open(targ_dir+'/OUTCAR','r') + tmp_mag = {} + it_mag = 0 + for line in datfile: + if 'NIONS = ' in line: + NIONS=int(line.split()[-1]) + if 'magnetization (x)' in line: + Null = datfile.readline() + Null = datfile.readline() + Null = datfile.readline() + tmp_list = [] + for i in range(int(NIONS)): + tmp_list.append(datfile.readline().split()[-1]) + tmp_mag[it_mag] = tmp_list + it_mag +=1 + datfile.close() + tot_mag = tmp_mag[it_mag-1] + mag_lines = '' + for i in range(len(tot_mag)): + mag_lines += str(' 1*'+tot_mag[i]) + return mag_lines diff --git a/src/module_subr.py b/src/module_subr.py index 436b286..8d74c6a 100644 --- a/src/module_subr.py +++ b/src/module_subr.py @@ -3,6 +3,7 @@ import math from itertools import product from itertools import combinations +from module_vector import dir_to_cart # This function is for calculate dot product of two vector def dotproduct(v1, v2): @@ -248,5 +249,59 @@ def set_atom(prim_at, M): # print sup_at_fin return sup_at_fin +def mkPairlist(atom_pos, axis, mag_true, cutoff): + list_pair = [] + mag_atom_pos = [] + for i in range(len(atom_pos)): + if mag_true[i] == 1: + mag_atom_pos.append(atom_pos[i]) + for i in range(len(mag_atom_pos)): + for j in range(len(mag_atom_pos)): + tmp_dist = calcDistwithPBC(axis, mag_atom_pos[i][0:3], mag_atom_pos[j][0:3], cutoff) + for dist in tmp_dist: + list_pair.append([i,j,str(dist)]) + return list_pair +def calcDistwithPBC(lparam, ptA, ptB, cutoff): + dist = [] + vecAB = np.array(ptA) - np.array(ptB) + for i in [-2,-1,0,1,2]: + for j in [-2,-1,0,1,2]: + for k in [-2,-1,0,1,2]: + vec_cur = np.copy(vecAB) + vec_cur = vec_cur + np.dot(np.array([i,j,k]), np.array(lparam)) + dist_cur = round(np.linalg.norm(vec_cur), 2) + + if dist_cur < cutoff and dist_cur > 0.1: + dist.append(dist_cur) + + return dist + +def check_pair(basis_mat,at_prim,mag_atom_list): + tmp = [] + for i in range(len(basis_mat)): + at_tmp = set_atom(at_prim, basis_mat[i][0]) + axis = basis_mat[i][1] + atom_pos_cart = [] + for j in range(len(at_tmp)): + atom_pos_cart.append(dir_to_cart(at_tmp[j][0:3],axis)+at_tmp[j][3:]) + mag_true = [] + mag_sum = 0 + for j in range(len(at_tmp)): + if at_tmp[j][3] in mag_atom_list: + mag_true.append(1) + mag_sum = mag_sum+1 + else: + mag_true.append(0) + info_pair = mkPairlist(atom_pos_cart,axis,mag_true, 5.0) + overlap_pair = [] + for k in range(len(info_pair)): + for j in range(len(info_pair)): + a,b,l_a = info_pair[k][1],info_pair[k][0],info_pair[k][2] + c,d,l_c = info_pair[j][1],info_pair[j][0],info_pair[j][2] + if int(a) == int(c) and int(b) == int(d) and float(l_a) != float(l_c): + overlap_pair.append(info_pair[k]) + if len(overlap_pair) == 0 : + tmp.append(basis_mat[i]) + return tmp diff --git a/src/module_vasprun.py b/src/module_vasprun.py index 1315af0..f87b9a2 100644 --- a/src/module_vasprun.py +++ b/src/module_vasprun.py @@ -156,7 +156,7 @@ def electronic_step_convergence_check(target): if elec_step == int(pygrep('NELM',target+'/OUTCAR',0,0).split(';')[0].split()[2]): make_amp2_log(target,'Electronic step is not converged.') algo = pygrep('ALGO',target+'/INCAR',0,0).split()[2] - if algo == 'Normal' or algo == 'All': + if algo == 'Normal': make_amp2_log(target,'Current ALGO is '+algo+' but it is not converged.') if spin == '2': if pygrep('BMIX_MAG',target+'/OUTCAR',0,0).split()[-1] == '1.00': @@ -168,6 +168,17 @@ def electronic_step_convergence_check(target): return 2 else: return 2 + if algo == 'All': + make_amp2_log(target,'Current ALGO is '+algo+' but it is not converged.') + if spin == '2': + if "F" in pygrep('LCHARG',target+'/OUTCAR',0,0).split("=")[1].split()[0].upper(): + make_amp2_log(target,'CHGCAR is saved to be used as a better initial guess') + wincar(target+'/INCAR',target+'/INCAR',[['LCHARG','T']],[]) + return 1 + else: + return 2 + else: + return 2 elif algo == 'Damped': wincar(target+'/INCAR',target+'/INCAR',[['ALGO','All']],[]) make_amp2_log(target,'ALGO changes from '+algo+' to All.') @@ -195,7 +206,22 @@ def electronic_step_convergence_check_CHGCAR(target): return 1 else: return 0 - +# check electronic step convergence for CHGCAR_conv +def electronic_step_convergence_check_CHGCAR_conv(target): + # 0: well converged, 1: can try to read CHGCAR + with open(target+'/OSZICAR','r') as inp: + fr_log = inp.readlines()[1:] + spin = pygrep('ISPIN',target+'/OUTCAR',0,0).split()[2] + elec_step = 0 + for ll in fr_log: + if ll.split()[0] == '1': + break + elec_step = elec_step+1 + if elec_step == int(pygrep('NELM',target+'/OUTCAR',0,0).split(';')[0].split()[2]): + wincar(target+'/INCAR',target+'/INCAR',[['NSW','0'],['NELM','100'],['LCHARG','T'],['ICHARG','1']],[]) + return 1 + else: + return 0 # check magnetism from relaxation def check_magnet(dir_in,min_mom): nion = int(pygrep('NION',dir_in+'/OUTCAR',0,0).split()[-1]) diff --git a/src/relax.py b/src/relax.py index ff94148..2c770b9 100644 --- a/src/relax.py +++ b/src/relax.py @@ -1,6 +1,6 @@ ########################################### -### Date: 2018-12-06 ### -### yybbyb@snu.ac.kr ### +### Date: 2020-11-05 ### +### mk01071@snu.ac.kr ### ########################################### import shutil, os, sys, subprocess, yaml from module_log import * @@ -72,7 +72,11 @@ make_amp2_log(dir_relax,'New calculation') copy_input(dir+'/INPUT0',dir_relax,POT) nsw = set_nsw(dir_relax+'/POSCAR',dir_relax+'/INCAR') - wincar(dir_relax+'/INCAR',dir_relax+'/INCAR',[['LCHARG','.F.']],[]) + if os.path.isfile(dir+'/INPUT0/CHGCAR_conv'): + shutil.copyfile(dir+'/INPUT0/CHGCAR_conv',dir_relax+'/CHGCAR') + wincar(dir_relax+'/INCAR',dir_relax+'/INCAR',[['LCHARG','.T.'],['ICHARG','1']],[]) + else: + wincar(dir_relax+'/INCAR',dir_relax+'/INCAR',[['LCHARG','.F.']],[]) converge_condition = set_ediffg(dir_relax+'/POSCAR',inp_rlx['force'],inp_rlx['pressure'],inp_rlx['energy']) if not converge_condition == 0: wincar(dir_relax+'/INCAR',dir_relax+'/INCAR',[['EDIFFG',str(converge_condition)]],[]) diff --git a/src/rerun_for_metal.py b/src/rerun_for_metal.py index 0b94054..b077f1a 100644 --- a/src/rerun_for_metal.py +++ b/src/rerun_for_metal.py @@ -11,7 +11,7 @@ from module_band import check_half_metal from input_conf import set_on_off from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2020-05-12' +code_data = 'Version '+__version__+'. Modified at 2020-01-15' # input from shell inp_file = sys.argv[1] @@ -28,24 +28,11 @@ Done_path = inp_yaml['directory']['done'] # Check metal in HSE -for pot_type in inp_yaml['hybrid_oneshot']['potential_type']: - if isinstance(pot_type,list): - if len(pot_type) == 1: - pot_cell = pot_type[0] - pot_point = pot_type[0] - else: - pot_cell = pot_type[0] - pot_point = pot_type[1] - else: - pot_cell = pot_type - pot_point = pot_type - if os.path.isfile(target+'/hybrid'+pot_cell+'_'+pot_point+'/Band_gap.log'): - with open(target+'/hybrid'+pot_cell+'_'+pot_point+'/Band_gap.log','r') as inp: - gap_log = inp.readline() - if not 'etal' in gap_log: - sys.exit() +if os.path.isfile(target+'/HSE/Band_gap.log'): + with open(target+'/HSE/Band_gap.log','r') as inp: + gap_log = inp.readline() # Check metal in GGA or LDA -if os.path.isfile(target+'/band_GGA/Band_gap.log'): +elif os.path.isfile(target+'/band_GGA/Band_gap.log'): with open(target+'/band_GGA/Band_gap.log','r') as inp: gap_log = inp.readline() elif os.path.isfile(target+'/band_LDA/Band_gap.log'): @@ -100,6 +87,9 @@ # Modify INCAR make_incar(target+'/INPUT0/POSCAR',target,src_path,inp_yaml['cif2vasp']['max_nelm']) +if os.path.isfile(target+'/INPUT0/CHGCAR_conv'): #Calculation with U would not have convergence problem + os.remove(target+'/INPUT0/CHGCAR_conv') + calc_out = 0 if set_on_off(cal_dic['kp_test']) == 1: From aad270aad6fa5714dbf233e58d17636a24ce771a Mon Sep 17 00:00:00 2001 From: Changho Hong Date: Tue, 1 Jun 2021 11:10:15 +0900 Subject: [PATCH 2/3] File trimmed --- src/config.yaml | 164 ------------------------------------------------ 1 file changed, 164 deletions(-) delete mode 100755 src/config.yaml diff --git a/src/config.yaml b/src/config.yaml deleted file mode 100755 index 3794516..0000000 --- a/src/config.yaml +++ /dev/null @@ -1,164 +0,0 @@ -directory: - submit: ./Submit - output: ./Output - done: ./Done - error: ./ERROR - src_path: ./src - pot_path_gga: /data/vasp4us/pot/PBE52 - pot_path_lda: /data/vasp4us/pot/LDA52 - -program: - vasp_std: ~/Initial/vasp/vasp_std - vasp_gam: /data/vasp4us/vasp_gam - vasp_ncl: /data/vasp4us/vasp_ncl - gnuplot: /gnuplot - mpi_command: mpirun - -calculation: - kp_test: 1 - encut_test: 1 - relaxation: 1 - band: 0 - density_of_states: 0 - hse_oneshot: 0 - dielectric: 0 - effective_mass: 0 - plot: 1 - magnetic_ordering: 1 - large_off: -1 - -cif2vasp: - pot_name: - GGA: - Ce: Ce_3 - LDA: -# Ge: Ge - soc_target: -# - Bi -# - Tl -# - Pb - u_value: -# La: 7.5 -# Ce: 4.5 -# All: 0.0 - magmom: 5 - max_nelm: 100 - incar: - -vasp_parallel: - npar: 2 - kpar: 1 - -convergence_test: - enconv: 0.01 - prconv: 10 - foconv: -1 - initial_kpl: 1 - max_kpl: 20 - enstart: 200 - enstep: 50 - enmax: 1000 - potential_type: GGA # Choose on potential in GGA, LDA, HSE - incar: - -relaxation: - potential_type: - - GGA -# - LDA -# - HSE - max_iteration: 10 - converged_ionic_step: 5 - pos_warning_percent: 1 - energy: -1 - pressure: 10 - force: 0.02 -# soc_on: 1 - incar: -# EDIFF: 1e-08 - -band_calculation: - kspacing_for_band: 0.005 - type_of_kpt: all # band (only for band structure plot) - y_min: 3 # set plot region, downward from VBM - y_max: 2 # set plot region, upward from CBM - potential_type: - - GGA -# - LDA -# - HSE - relax_check: 1 - incar: - KPAR: 2 - NPAR: 2 - -hybrid_oneshot: - alpha: 0.25 - potential_type: # The potential used for lattice parameter optimization - - GGA # If one variable is inserted, we use the lattice parameter and the points of VBM and CBM with that potential. -# - - HSE # If two variables are inserted, we use the lattice parameter with above potential and the points of VBM and CBM with below potential. -# - GGA - fermi_width: 0.3 - vb_dos_min: 1 - vb_dos_max: 3 - energy_width_for_extreme: 0.5 - search_space_for_extreme: 0.05 - cutoff_df_dvb: 0.3 - band_structure_correction: 1 - incar: - -density_of_states: - potential_type: - - GGA -# - LDA -# - HSE - relax_check: 1 - kp_multiplier: 2 - y_min: 3 # set plot region, downward from VBM - y_max: 2 # set plot region, upward from CBM - incar: - -dielectric: - kp_multiplier: 2 - potential_type: - - GGA -# - LDA - incar: - EDIFF: 1e-08 -# EDIFFG: 1e-04 - metal_check: 1 - relax_check: 1 - -effective_mass: - carrier_type: - - hole - - electron - pocket_search_dist: 0.02 - grid_size_for_searching: 0.02 - grid_size_for_calculation: 0.02 - temperature_for_fermi: 300 - fermi_for_cutoff: 0.99 - potential_type: - - GGA -# - LDA -# - HSE - incar: - -magnetic_ordering: - potential_type: GGA # Choose on potential in GGA, LDA, HSE (If from relax:0, it must be same to potential_type in convergence.) - minimum_moment: 0.5 - from_relax: 1 - cutoff_for_parameter: 5.0 - tolerance: 0.01 # Tolerance to define identical pair for Ising coefficient - energy_tolerance: 0.00001 # Tolerance to find lowest energy spin ordering candidates - maximum_atoms_number: -1 - maximum_pair_type: 20 - genetic_algorithm: - min_iteration: 0 - max_iteration: -1 - selection_coeff: 2 - population: - best: 15 - crossover: 30 - mutation: 30 - random: 15 - incar: -# NPAR: 2 From d7fd33568819988587a670bfedf166c89ed41119 Mon Sep 17 00:00:00 2001 From: Changho Hong Date: Tue, 1 Jun 2021 11:13:09 +0900 Subject: [PATCH 3/3] File trimmed 2 --- src/band.py | 9 +++++++-- src/cutoff.py | 9 +++++++-- src/dos.py | 9 +++++++-- src/hse_gap.py | 23 +++++++++++++++++------ src/kpoint.py | 9 +++++++-- src/magnetic_ordering.py | 2 +- src/rerun_for_metal.py | 21 +++++++++++++++++---- 7 files changed, 63 insertions(+), 19 deletions(-) diff --git a/src/band.py b/src/band.py index 9fe0ffd..1e5ff31 100644 --- a/src/band.py +++ b/src/band.py @@ -10,7 +10,7 @@ from module_hse import * from input_conf import set_on_off from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2019-12-17' +code_data = 'Version '+__version__+'. Modified at 2020-11-05' # Set input dir = sys.argv[1] @@ -209,9 +209,14 @@ # Drawing band structure (nkpt_for_band is used for hse band structure.) plot_band_structure(spin,[x[-1*nkpt_for_band:] for x in Band],fermi,dir_band+'/xtic.dat',dir_band+'/xlabel.dat',[inp_band['y_min'],inp_band['y_max']],dir_band) +if not os.path.isfile(inp_yaml['program']['gnuplot']): + make_amp2_log(dir_band,'If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir_band) - subprocess.call([gnuplot,dir_band+'/band.in']) + try: + subprocess.call([gnuplot,dir_band+'/band.in']) + except: + make_amp2_log(dir_band,'Error occured drawing figure. Please check gnuplot.') with open(dir_band+'/amp2.log','r') as amp2_log: with open(dir+'/amp2.log','a') as amp2_log_tot: diff --git a/src/cutoff.py b/src/cutoff.py index 5e4763a..0b9304b 100644 --- a/src/cutoff.py +++ b/src/cutoff.py @@ -10,7 +10,7 @@ from module_converge import * import math from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2019-12-17' +code_data = 'Version '+__version__+'. Modified at 2020-11-05' # Set input dir = sys.argv[1] @@ -201,9 +201,14 @@ wincar(dir+'/INPUT0/INCAR',dir+'/INPUT0/INCAR',[['ENCUT',str(ENCUT-3*ENSTEP)]],[]) make_conv_dat(dir+'/cutoff','cutoff') +if not os.path.isfile(inp_yaml['program']['gnuplot']): + make_amp2_log(dir+'/cutoff','If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir+'/cutoff') - subprocess.call([gnuplot,dir+'/cutoff/conv_plot.in']) + try: + subprocess.call([gnuplot,dir+'/cutoff/conv_plot.in']) + except: + make_amp2_log(dir+'/cutoff','Error occured drawing figure. Please check gnuplot.') with open(dir+'/cutoff/amp2.log','r') as amp2_log: with open(dir+'/amp2.log','a') as amp2_log_tot: diff --git a/src/dos.py b/src/dos.py index 8d70f31..094d41e 100644 --- a/src/dos.py +++ b/src/dos.py @@ -9,7 +9,7 @@ from module_dos import * from input_conf import set_on_off from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2019-12-17' +code_data = 'Version '+__version__+'. Modified at 2020-11-05' # Set input dir = sys.argv[1] @@ -247,9 +247,14 @@ make_dos_in(dir_dos,atom_name,spin,len(par_dos[0][0]),[inp_dos['y_min'],inp_dos['y_max']+gap]) +if not os.path.isfile(inp_yaml['program']['gnuplot']): + make_amp2_log(dir_dos,'If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir_dos+'/Pdos_dat') - subprocess.call([gnuplot,dir_dos+'/Pdos_dat/dos.in']) + try: + subprocess.call([gnuplot,dir_dos+'/Pdos_dat/dos.in']) + except: + make_amp2_log(dir_dos,'Error occured drawing figure. Please check gnuplot.') make_amp2_log(dir_dos,'DOS calculation is done.') diff --git a/src/hse_gap.py b/src/hse_gap.py index 71a209c..df40077 100644 --- a/src/hse_gap.py +++ b/src/hse_gap.py @@ -1,6 +1,7 @@ ########################################### -### Date: 2018-12-05 ### +### Date: 2020-11-05 ### ### yybbyb@snu.ac.kr ### +### mk01071@snu.ac.kr ### ########################################### # This is for estimating band gap with PBE@HSE scheme. import shutil, os, sys, subprocess, yaml @@ -12,7 +13,7 @@ from input_conf import set_on_off from module_relax import * from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2019-12-17' +code_data = 'Version '+__version__+'. Modified at 2020-11-05' # Set input dir = sys.argv[1] @@ -171,9 +172,9 @@ incar_for_hse(dir_hse+'/INCAR') if PBE0_on == 1: inp_hse['alpha'] = calc_alpha_auto(diel_path+'/dielectric.log') - wincar(dir_hse+'/INCAR',dir_hse+'/INCAR',[['NSW','0'],['HFSCREEN','0.0'],['ALGO','ALL'],['AEXX',str(inp_hse['alpha'])]],[]) + wincar(dir_hse+'/INCAR',dir_hse+'/INCAR',[['NSW','0'],['HFSCREEN','0.0'],['ALGO','All'],['AEXX',str(inp_hse['alpha'])]],[]) else: - wincar(dir_hse+'/INCAR',dir_hse+'/INCAR',[['NSW','0'],['HFSCREEN','0.2'],['ALGO','ALL'],['AEXX',str(inp_hse['alpha'])]],[]) + wincar(dir_hse+'/INCAR',dir_hse+'/INCAR',[['NSW','0'],['HFSCREEN','0.2'],['ALGO','All'],['AEXX',str(inp_hse['alpha'])]],[]) incar_from_yaml(dir_hse,inp_hse['incar']) mag_on = check_magnet(dir+'/relax_'+pot_cell,inp_yaml['magnetic_ordering']['minimum_moment']) vasprun = make_incar_for_ncl(dir_hse,mag_on,kpar,npar,vasp_std,vasp_gam,vasp_ncl) @@ -219,9 +220,14 @@ for k in range(len(KPT)): Band[n][k][i] = Band[n][k][i] + E_shift plot_band_corrected_structure(spin,Band,eVBM,dir_band+'/xtic.dat',dir_band+'/xlabel.dat',[inp_band['y_min'],inp_band['y_max']+E_shift],dir_band) + if not os.path.isfile(inp_yaml['program']['gnuplot']): + make_amp2_log(dir_hse,'If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir_band) - subprocess.call([gnuplot,dir_band+'/band_corrected.in']) + try: + subprocess.call([gnuplot,dir_band+'/band_corrected.in']) + except: + make_amp2_log(dir_hse,'Error occured drawing figure. Please check gnuplot.') # Metal in GGA. Band reordering is required. else: @@ -252,9 +258,14 @@ fermi = get_fermi_level(Band_reorder,nelect,ncl) if calc_gap(fermi,spin,ncl,KPT[:num_kpt_for_image],Band_reorder,nelect) > 0.01: plot_band_corrected_structure(spin,Band_reorder,eVBM,dir_band+'/xtic.dat',dir_band+'/xlabel.dat',[inp_band['y_min'],inp_band['y_max']+float(gap)],dir_band) + if not os.path.isfile(inp_yaml['program']['gnuplot']): + make_amp2_log(dir_hse,'If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir_band) - subprocess.call([gnuplot,dir_band+'/band_corrected.in']) + try: + subprocess.call([gnuplot,dir_band+'/band_corrected.in']) + except: + make_amp2_log(dir_hse,'Error occured drawing figure. Please check gnuplot.') else: make_amp2_log(dir_hse,'Warning. We cannot correct the band structure.') else: diff --git a/src/kpoint.py b/src/kpoint.py index 88bc8e9..cc28a3b 100644 --- a/src/kpoint.py +++ b/src/kpoint.py @@ -9,7 +9,7 @@ from module_vasprun import * from module_converge import * from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2019-12-17' +code_data = 'Version '+__version__+'. Modified at 2020-11-05' # Set input dir = sys.argv[1] @@ -222,9 +222,14 @@ wincar(dir+'/INPUT0/INCAR',dir+'/INPUT0/INCAR',[["ALGO","All"]],[]) make_conv_dat(dir+'/kptest','kpoint') +if not os.path.isfile(inp_yaml['program']['gnuplot']): + make_amp2_log(dir+'/kptest','If you want to draw figure, please check the path of gnuplot.') if inp_yaml['calculation']['plot'] == 1: os.chdir(dir+'/kptest') - subprocess.call([gnuplot,dir+'/kptest/conv_plot.in']) + try: + subprocess.call([gnuplot,dir+'/kptest/conv_plot.in']) + except: + make_amp2_log(dir+'/kptest','Error occured drawing figure. Please check gnuplot.') with open(dir+'/kptest/amp2.log','r') as amp2_log: with open(dir+'/amp2.log','a') as amp2_log_tot: diff --git a/src/magnetic_ordering.py b/src/magnetic_ordering.py index b0dc5c8..c697fda 100644 --- a/src/magnetic_ordering.py +++ b/src/magnetic_ordering.py @@ -13,7 +13,7 @@ from module_relax import * from input_conf import set_on_off from _version import __version__ -code_data = 'Version '+__version__+'. Modified at 2019-12-17' +code_data = 'Version '+__version__+'. Modified at 2020-11-05' # Set input dir = sys.argv[1] diff --git a/src/rerun_for_metal.py b/src/rerun_for_metal.py index b077f1a..c368ce9 100644 --- a/src/rerun_for_metal.py +++ b/src/rerun_for_metal.py @@ -28,11 +28,24 @@ Done_path = inp_yaml['directory']['done'] # Check metal in HSE -if os.path.isfile(target+'/HSE/Band_gap.log'): - with open(target+'/HSE/Band_gap.log','r') as inp: - gap_log = inp.readline() +for pot_type in inp_yaml['hybrid_oneshot']['potential_type']: + if isinstance(pot_type,list): + if len(pot_type) == 1: + pot_cell = pot_type[0] + pot_point = pot_type[0] + else: + pot_cell = pot_type[0] + pot_point = pot_type[1] + else: + pot_cell = pot_type + pot_point = pot_type + if os.path.isfile(target+'/hybrid'+pot_cell+'_'+pot_point+'/Band_gap.log'): + with open(target+'/hybrid'+pot_cell+'_'+pot_point+'/Band_gap.log','r') as inp: + gap_log = inp.readline() + if not 'etal' in gap_log: + sys.exit() # Check metal in GGA or LDA -elif os.path.isfile(target+'/band_GGA/Band_gap.log'): +if os.path.isfile(target+'/band_GGA/Band_gap.log'): with open(target+'/band_GGA/Band_gap.log','r') as inp: gap_log = inp.readline() elif os.path.isfile(target+'/band_LDA/Band_gap.log'):