class0_incar.py

import numpy as np 

''' 
 Given: directory, a dictionary of input values 
 Output: INCAR in the given directory 
''' 

def write_INCAR(direc, dictionary={}, incarname='INCAR'): 
	if direc[-1]!='/':
		direc+='/'
	comments = {'SYSTEM':'# Title\n\n#  (Re)start parameters:\n', 
			'ISTART':'# read WAVECAR?: 0-no begin from scratch, 1-yes restart with ionstant energy cutoff and Orbitals are read from WAVECAR, 2-yes constant basis, Use 0 (if no WAVECAR) or 1 (If have WAVECAR to save time)\n', 
			'ICHARG':'# how VASP constructs the initial charge density. Read CHGCAR?: 2-no if ISTART=0 (superpositions of atomic charge densities), 0-no (Calculate charge density from initial wave function in WAVECAR), 1-yes,read and extropolate from old positions (on CHGCAR) to new positions using linear combination of charge densities, self-consistent calculations, 11-yes, non-self-consisitent calculations. Charge density will be kept constant during the whole electronic minimization. To obtain the eigenvalues (for band structure plots) or the DOS for a given charge density read from CHGCAR. The selfconsistent CHGCAR file must be determined beforehand doing by a fully self-consistent calculation with a k-point grid spanning the entire Brillouin zone.\n\n#  Parallelization\n', 
			'KPAR':'# number of k-points that are to be treated in parallel\n', 
			'NCORE':'# number of compute cores that work on an individual orbital: NCORE =number-of-cores / NPAR\n', 
			'NPAR':'# number of nodes for IALGO=48, otherwise 1\n', 
			'LPLANE':'# by plane data distribution (FFT), set to true if neywork is slow\n', 
			'LSCALU':'# whether LU decomposition is done in parallel\n\n#  Accuracy\n', 
			'EDIFF':'# convergence condition\n', 
			'EDIFFG':'# negative: force convergence condition\n', 
			'ENCUT':'# specifies the cutoff energy for the planewave basis set in eV.\n', 
			'PREC':'# Low | Normal | Single | Accurate \n', 
			'ADDGRID':'# whether grid for augmentation charges should have 8 times more grid points that NG(X,Y,x)F grid. For accurate calcualtions use TRUE\n', 
			'LMAXMIX':'# l quantum number up to whcih onsite PAW charge densities are passsed through the charge density mixer and written to WAVECAR (for LDA+U use 4 for d systems and 6 for f systems)\n\n# Electronic Steps\n', 
			'NELECT':'# total number of electrons\n', 
			'NELM':'# max number of electronic steps\n', 
			'NELMIN':'# forces min number (=NELMIN) of electronic steps (increase for relaxation: 5 for inaccurate relaxation, 8 for accurate) to guarantee that the forces are well converged at each step\n', 
			'NELMDL':'# number of eectronic steps at the beginning of the first ionic relazation for which charge density is not upated (absolute value of course) \n\n#  Ionic steps\n', 
			'NSW':'# number of ionic steps\n', 
			'ISIF':'# determines whether the stress tensor is calculated and which degrees-of-freedom (ionic positions, cell volume, and cell shape) are allowed to change.. what to relax?, find stress tensor?: 0 (for IBRION=0 molecular dynamics) only ions,no, 1 only ions can move, only trace, 2 only ions, fix shape and volume, used for defect calculation, yes, 3 ions,shape,vol, used in relaxation to find ground state energy, yes, 4 ions,shape, yes, 5 shape,yes, 6 shape,vol,yes, 7 vol can, not shape nor, yes\n', 
			'IBRION':'# how the ions are updated and moved: -1 ions fixed (NSW>0 does not make much sense), 0-molecular dynamics Verlet algorithm (only constant volume MD\'s are possible), 1 ion relaxation (quasi-Newton algorithm, relax the ions into their instantaneous groundstate, fails badly if the initial positions are a bad guess, use IBRION=2 for this case), 2 ion relaxation (conjugate gradient algorithm), 3 ion relaxation (damped molecular dynamics, damping factor set by SMASS), 5 Gamma point phonons\n', 
			'POTIM':'# scalling factor for forces if relaxation is done, time step (fs) for MD,for IBRION=5 step size for displacement of ions in froze phonon calcualtions\n', 
			'IWAVPR':'# predition, how orbitals/densities are extrapolated from one ionic configuration to the next configuration: 0-static no prediction, 1-charge, 2-wave, 3-Comb\n', 
			'LCORR':'# calculate Harris-correction to forces\n\n# BZ integration\n', 
			'ISMEAR':'# how partial occupancies are found: N > 0 use the method N-order Methfessel-Paxton; 0 Gaussian(sigma), -1 Fermi(sigma), -2 partial occupancies are read in from the WAVECAR or INCAR file, and kept fixed throughout run, -5 tetrahedron with corrections\n', 
			'SIGMA':'# broadening width [eV] \n\n\n# HSE setting  \n', 
			'KSPACING':'#\n', 
			'GGA':'#\n', 
			'LHFCALC':'# HF type calculations are performed \n', 
			'HFSCREEN':'# truncate the long range Fock potential\n', 
			'ALGO':'# Default:Normal. specify the electronic minimisation algorithm.  Damped|All|Normal|Exact|Fast  !Normal uses blocked Davidson iteration scheme !Exact for Exact Diagonalization !Damped for metal & small-gap semiconductors !All selects an all band simultaneous update of orbitals !Fast selects a faily robust mixture of the Davidson and RMM-DIIS algorithms.\n', 
			'TIME':'#\n', 
			'NKPED':'#\n', 
			'PRECFOCK':'#\n', 
			'AEXX':'#\n', 
			'AGGAX':'#\n', 
			'AGGAC':'#\n', 
			'ALDAC':'#\n\n# wannier output\n', 
			'LWANNIER90':'# F is default. T to have wannier input files\n\n\n# Zero field splitting\n', 
			'LDMATRIX':'# Will give D matrix: Spin-spin contribution to zero-field splitting tensor (MHz). not documented anywhere\n\n# Hyperfine parameter\n', 
			'LHYPERFINE':'#  compute the hyperfine tensors at the atomic sites. The total Fermi contact term is given by Atot. Beware: for the moment we have chosen NOT to include the core contributions A1c in Atot. If you so want, you should add them by hand to Atot.\n', 
			'NGYROMAG':'# specify nuclear gyromagnetic ratios. Need to read which nuclear isotope is used in experiment, use the gamma for that isotope. Specify one number for each of the N species on the POSCAR file. If not set NGYROMAG in the INCAR file, VASP assumes a factor of 1 for each species.\n\n#  Magnetic (Non-collinear) calculaton:\n', 
			'ISPIN':'# Magnetism: 1 nonmagnetic, 2 magnetic and PROCAR contains a second set of projections for the spin down channel.\n', 
			'MAGMOM':'# Initial magnetic moment ( for noncollinear three componens of magnetic moment)for each atom (order as in POSCAR);\n', 
			'VOSKOWN':'# Vosko\'s interpolation for the exchange correlation potential\n', 
			'NUPDOWN':'# -1:not set   ;;  fixed spin moment\n', 
			'LNONCOLLINEAR':'# noncollinear magnetism\n', 
			'LSORBIT':'# spin-orbit coupling\n', 
			'SAXIS':'# spin quantization axis\n', 
			'LORBMOM':'# whether the orbital moments are written out or not. must have LSORBIT=T\n\n# dielectric constant calculation: # PEAD (Perturbation Expression After Discretization): determine the ground state of an insulating system under the application of a finite homogeneous electric field\n', 
			'LPEAD':'#As an alternative to solving a linear Sternheimer equation (one method to calculate displacements, I think), one may compute the derivative of the cell-periodic part of the orbitals w.r.t. k  from finite differences by specifying LPEAD=.TRUE.\n', 
			'LEPSILON':'# determines the static matrix, does not support HSE\n', 
			'LOPTICS':'# computes frequency dependent dielectric matrix after the ground state has been determined\n', 
			'CSHIFT':'# CSHIFT=0.1 by default in linear response calculations, and otherwise in GW. CSHIFT sets the (small) complex shift η in the Kramers-Kronig transformation. If the gap is very small (i.e. approaching two times CSHIFT), slight inaccuracies in the static dielectric constant are possible, which can be remedied by decreasing CSHIFT. If CSHIFT is further decreased, it is strongly recommended to increase the parameter NEDOS to values around 2000.\n', 
			'LCALCEPS':'# the macroscopic ion-clamped static dielectric tensor, Born effective charge tensors, and the ion-clamped piezoelectric tensor of the system are determined from the reponse to finite electric fields\n\n# Constrained occupation calculation: #set ISMEAR:-2 \n', 
			'FERWE':'# FERWE sets the occupancies of the states in the up-spin channel.the left number is the number of states to be filled with the right number which is the occupation, so in both spin up and spin down, we fill the first 300 bands with 1 electron (the valence band).. in spin up, we then fill the next 3 states with 1 electron each, in spin down, we only put one electron in the middle state (corresponding to lifting the one down electron up one level).. then we have 999 empty bands (the conduction band) #300 is probably not the right number, you should figure this out from your EIGENVAL file, 999 is probably also not right, but VASP does not care if there are more empty bands than are necessary\n', 
			'FERDO':'# FERDO sets the occupancies of the states in the down-spin channel for ISMEAR:-2 and ISPIN:2. Start by using ALGO: Normal, if this fails, try ALGO: All, if this fails, try ALGO: Damped,\n\n# Electronic part\n', 
			'IALGO':'# Hamiltonian diagonalisation algorithm: 8 conjugate gradient, 38 Davidson scheme, 48 residual minimization\n', 
			'NSIM':'# number of bands that are optimized at the same time for IALGO=48\n', 
			'LDIAG':'#sub-space diagonalisation\n', 
			'NBANDS':'# number of bands\n', 
			'LREAL':'# whether scalar product of projectros and the wavefunctions is evaluated in real space (as opposed to k-space) (.F. for small systems; AUTO for large).\n\n# VdW correction\n', 
			'LVDW':'# turn on VdW correction\n', 
			'IVDW':'# 11 and 12 for vdW D3 correction\n\n# Symmetry\n', 
			'ISYM':'# symmetry: -1 the use of symmetry is switched off completely, 0-NO,reduces the sampling of the Brillouin zone compared to -1, should be set for molecular dynamics, i.e. IBRION=0., 1-yes, 2-yes efficiently\n', 
			'SYMPREC':'# 1E-5 is default. Determines how accurately the positions in the POSCAR file must be specified\n\n# Pressure\n', 
			'PSTRESS':'# The unit of PSTRESS is kB (1GPa=10kB). During ionic minimization, an energy term E=V*PSTRESS} is added to the total energy and the pressure set by PSTRESS is added to the diagonals of the stress tensor. This allows performing structure optimization at a specific external pressure.\n\n# DOS calcualtions: DOSCAR\n', 
			'NEDOS':'# number of energies for DOS\n', 
			'EMIN':'# min energy [eV] for DOS\n', 
			'EMAX':'# max energy [eV] for DOS\n', 
			'RWIGS':'# Wigner-Seitz radii, one for each type of atom (order as in POSCAR), if provided partial DOS and l- and site- projected wavefunctions are calcualted.\n\n#  Which files to be written\n', 
			'LWAVE':'# print WAVECAR, True by default\n', 
			'LCHARG':'# print CHGCAR, True by default\n', 
			'LVTOT':'# print LOCPOT: plot both Hartree and exchange correlation. TOT means total.\n', 
			'LVHAR':'# print LOCPOT (electrostatic potential). The term HAR means only hartree fock term. Used for Freysoldt correction.\n', 
			'LORBIT':'# print PROCAR or PROOUT. Project charge density onto each atom. 11: project onto spherical harmonics, generate lm-dependent PROCAR; whether to produce lm (angular momentum quantum number) resolved DOSCAR, PORCAR, PROOUT files (see manual); # Set ISPIN:2 for separate spin projections. PROCAR contains a second set of projections for the spin down channel\n', 
			'LPARD':'# print PARCHG. True: evaluate partial charge densities. Orbitals read from WAVECAR must be converged in a separate prior run. Need NSW:0, IBRION:-1, ISTART:1, LPARD:.True., IBAND:457 \n', 
			'IBAND':'# For PARCHG. Specify the band for viewing partial charge densities. \n', 
			'LSEPB':'# For PARCHG. calculate each band separately\n\n# Electronic mixing\n', 
			'IMIX':'# Mixing scheme: 1-Kerker, 2-Tchebycheff, 4-Broyden or Pulay (depending on WC)\n', 
			'WC':'# weight factor for each iterative step (only for IMIX=4): >0-Pullay; 0-Broyden\n', 
			'INIMIX':'# Type of initial mixing (only for IMIX=4): onic mixing, 0-linear(AMIX); 1-Kerker(AMIX,BMIX)\n', 
			'AMIN':'# Minimal mixing parameter\n', 
			'AMIX':'# Linear mixing parameter \n', 
			'BMIX':'# cutoff wavevector for Kerker method\n', 
			'AMIX_MAG':'# linear mixing paramter for mixing magneitzation density\n', 
			'BMIX_MAG':'# Kerker wavevector cutoff fot magnetization denisty mixing\n\n# LDA+U calculation:\n', 
			'LDAU':'# flag\n', 
			'LDAUTYPE':'# LDA+U type: 1 Liechtenstein,  2-spherically averaged (Dudarev)\n', 
			'LDAUL':'# orbitals: -1 no, 0 s, 1 p, 2 d, 3 f(vasp.4.6 only) for each species\n', 
			'LDAUU':'# values of U for each species\n', 
			'LDAUJ':'# values of J for each species\n', 
			'LDAUPRINT':'# \n', 
			'LASPH':'# Logical ASPHerical: .TRUE.-nonspherical contributions from the gradient corrections inside the PAW spheres will be included\n\n# Berry phase calculation:\n', 
			'LBERRY':'# flag\n', 
			'IGPAR':'# direction of integration in the reciprocal unit cell (1-3)\n', 
			'NPPSTR':'# number of points in this direction\n', 
			'DIPOL':'# cell center\n\n# Elastic band method\n', 
			'IMAGES':'# number of images\n', 
			'SPRING':'# spring constant (ansolute value of course)r(AMIX); 1-Kerker(AMIX,BMIX)\n', 
			'LCLIMB':'# the climbing image algorithm, a more accurate method of finding the saddle points\n', 
			'ICHAIN':'# Indicates which method to run. NEB (ICHAIN=0) is the default\n', 
			'IOPT':'# optimizer: 0-default VASP optimizers specified from IBRION, 1-LBFGS limited-memory Broyden-Fletcher-Goldfarb-Shanno, 2-conjugate gradient, 3-Quick-Min, 4-steepest descent, 7-fast inertial relaxation engine\n\n# Molecular dynamics \n## need IBRION:0, POTIM, NSW, ISIF:2 for NVT; 3 for NPT, TEBEG, MDALGO \n', 
			'SMASS':'# -3: NVE micro cannonical ensemble. -2 the velocities are kept constant as POTIM * POSCAR_velocity, used for frozen phonons, or dimers. SMASS=0, a Nosé-mass corresponding to period of 40 time steps will be chosen. -1: annealing, velocities are scaled each NBLOCK step: T=TEBEG+(TEEND-TEBEG)×NSTEP/NSW, NBLOCK controls how often the kinetic energy is scaled.\n', 
			'TEBEG':'# starting temperature (in K) for an ab-initio molecular dynamics run (IBRION=0) and other routines\n', 
			'TEEND':'# the final temperature (in K), equals TEBEG by default, for molecular-dynamics run (IBRION=0; SMASS=−1)\n', 
			'NBLOCK':'# After NBLOCK ionic steps the pair-correlation function and the DOS are calculated and write CONTCAR file. Usually set NBLOCK=1. Only for many 1000 steps or machine learning, NBLOCK=10 or even 100, to avoid large XDATCAR-files.\n', 
			'KBLOCK':'# After NBLOCK*KBLOCK steps, write XDATCAR file\n', 
			'MDALGO':'# 0: standard molecular dynamics, Should have the same results as MDALGO=2. 2:Nose-Hoover thermostat, add the kinetic and potential terms of the fictitious coordinates. 3: Langevin thermostat, NVT needs ISIF=2 while NPT needs ISIF=3\n', 
			'MAXMIX':'# specifies the maximum number of steps stored in the Broyden mixer\n', 
			'LANGEVIN_GAMMA':'# friction coefficients for atomic degree of freedom, separate coefficient for each type of element\n', 
			'LANGEVIN_GAMMA_L':'# friction coefficient (in ps-1) for lattice degrees-of-freedom in case of Parrinello-Rahman dynamics\n\n# Machine learning force field\n', 
			'ML_LMLFF':'# might be true to have MLFF\n', 
			'ML_MODE':'# run: only preditions no training. train: on-the-fly training, if ML_AB file is present, will restart the training. refit: refitting from ML_AB file and number of local reference confs is fixed. select: new MLFF generated from ML_AB, and local reference confs are auto chosen\n', 
			'ML_IALGO_LINREG':'# 1: algorithm is Bayesian linear regression. 4: Singular value decomposition with Tikhonov regularization, Usable with ML_MODE = REFIT\n', 
			'ML_SIONS':'# sigma_ions of Gaussian functions used for broadening the atomic distributions of the radial descriptor\n', 
			'ML_MRB2':'# number of radial basis functions used to expand the angular descriptor\n'}
	default_vals = {'SYSTEM':'bulk	', 
			'ISTART':'1	', 
			'ICHARG':'1	', 
			'KPAR':'4	', 
			'NCORE':'4	', 
			'NPAR':'8	', 
			'LPLANE':'F	', 
			'LSCALU':'F	', 
			'EDIFF':'0.00001	', 
			'EDIFFG':'-0.01	', 
			'ENCUT':'600	', 
			'PREC':'Accurate	', 
			'ADDGRID':'T	', 
			'LMAXMIX':'4	', 
			'NELECT':'473	', 
			'NELM':'160	', 
			'NELMIN':'5	', 
			'NELMDL':'-5	', 
			'NSW':'200	', 
			'ISIF':'3	', 
			'IBRION':'2	', 
			'POTIM':'0.5	', 
			'IWAVPR':'1	', 
			'LCORR':'T	', 
			'ISMEAR':'0	', 
			'SIGMA':'0.05	', 
			'KSPACING':'0.30 ', 
			'GGA':'HSE06     ', 
			'LHFCALC':'.TRUE.	', 
			'HFSCREEN':'0.2	', 
			'ALGO':'Normal ', 
			'TIME':'0.4	', 
			'NKPED':'2	', 
			'PRECFOCK':'Fast	', 
			'AEXX':'0.32	', 
			'AGGAX':'0.675	', 
			'AGGAC':'1.0	', 
			'ALDAC':'1.0	', 
			'LWANNIER90':'T ', 
			'LDMATRIX':'True ', 
			'LHYPERFINE':' .TRUE. ', 
			'NGYROMAG':' gamma_1  gamma_2 ... gamma_N ', 
			'ISPIN':'1	', 
			'MAGMOM':'2	', 
			'VOSKOWN':'1	', 
			'NUPDOWN':'-1	', 
			'LNONCOLLINEAR':'T	', 
			'LSORBIT':'T	', 
			'SAXIS':'0 0 1	', 
			'LORBMOM':'T	', 
			'LPEAD':'.False.	', 
			'LEPSILON':'.TRUE.	', 
			'LOPTICS':'.TRUE.	', 
			'CSHIFT':'0.1	', 
			'LCALCEPS':'.TRUE.	', 
			'FERWE':'300*1.0 1*1.0 1*1.0 1*1.0 999*0.0	', 
			'FERDO':'300*1.0 1*0.0 1*1.0 1*0.0 999*0.0	', 
			'IALGO':'38	', 
			'NSIM':'4	', 
			'LDIAG':'T	', 
			'NBANDS':'160	', 
			'LREAL':'Auto	', 
			'LVDW':'.True.	', 
			'IVDW':'11 ', 
			'ISYM':'-1	', 
			'SYMPREC':'1E-3	', 
			'PSTRESS':'0 ', 
			'NEDOS':'1001	', 
			'EMIN':'-7.0	', 
			'EMAX':'12.0	', 
			'RWIGS':'[real array]	', 
			'LWAVE':'F	', 
			'LCHARG':'F	', 
			'LVTOT':'T	', 
			'LVHAR':'.TRUE.	', 
			'LORBIT':'11	', 
			'LPARD':'.TRUE.	', 
			'IBAND':'11 12	', 
			'LSEPB':'.TRUE.  ', 
			'IMIX':'4	', 
			'WC':'1000	', 
			'INIMIX':'1	', 
			'AMIN':'0.1	', 
			'AMIX':'0.4	', 
			'BMIX':'1	', 
			'AMIX_MAG':'1.6	', 
			'BMIX_MAG':'1.0	', 
			'LDAU':'T	', 
			'LDAUTYPE':'2	', 
			'LDAUL':'2 2 -1	', 
			'LDAUU':'1.73 4.0 0.0	', 
			'LDAUJ':'0.70 0.58 0.0	', 
			'LDAUPRINT':'2	', 
			'LASPH':'T	', 
			'LBERRY':'F	', 
			'IGPAR':'3	', 
			'NPPSTR':'6	', 
			'DIPOL':'0.25 0.25 0.25	', 
			'IMAGES':'1	', 
			'SPRING':'-5	', 
			'LCLIMB':'.TRUE.	', 
			'ICHAIN':'0	', 
			'IOPT':'3	', 
			'SMASS':'0   ', 
			'TEBEG':'0   ', 
			'TEEND':'TEBEG   ', 
			'NBLOCK':'1  ', 
			'KBLOCK':'50 ', 
			'MDALGO':'2  ', 
			'MAXMIX':'50 ', 
			'LANGEVIN_GAMMA':' 0 0 0    ', 
			'LANGEVIN_GAMMA_L':' 0  ', 
			'ML_LMLFF':' .TRUE. ', 
			'ML_MODE':' train   ', 
			'ML_IALGO_LINREG':' 1   ', 
			'ML_SIONS':' 0.5    ', 
			'ML_MRB2':' 8   '}
	with open(direc+incarname, 'w') as f:
		for key in default_vals.keys():
			if key in dictionary:
				f.write('%s%s=%s %s' % ( '', key, str(dictionary[key]), comments[key]) )
			else:
				f.write('%s%s=%s%s' % ( '# ', key, default_vals[key], comments[key]) ) 
#write_INCAR('./')