-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsimulationparameters.py
179 lines (140 loc) · 7.07 KB
/
simulationparameters.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
#
# Routines for translating simulation parameters and k-point meshes
#
import numpy as np
from scipy.constants import physical_constants as constants
from pseudopotentials import qe_pppath as ppath
RytoeV = 1./constants["Rydberg constant times hc in eV"][0]
def set_espresso_param(parameters, kw, kwarg):
#
# Update espresso parameter dictionary with keyword and keyword argument
#
parameters[kw] = kwarg
return()
def vasp_to_espresso(vasp_params, basic=False, debug=True, version=6.7):
#
# Go through the vasp parameters and create dictionary of espresso parameters
#
# Initial (minimal set of params) to run simulation
espresso_params = {
# CONTROL
"calculation" : "scf",
"prefix" : "espresso",
"outdir" : "tmp",
"pseudo_dir" : ppath,
"verbosity" : "high",
"title" : "name",
"tprnfor" : True,
"disk_io" : "low",
# SYSTEM
"ecutwfc" : 40,
# ELECTRONS
}
if not basic:
# Basic option ignores all following vasp parameters, can be used to check symmetry etc.
for param in vasp_params:
if param == 'SYSTEM':
set_espresso_param(espresso_params, "title", vasp_params[param])
elif param == 'ALGO':
# Algorithm for solution of the SCF cycle
if (vasp_params[param] == "Normal"):
set_espresso_param(espresso_params, "diagonalization", "david")
elif (vasp_params[param] == "Fast") or (vasp_params[param] == "VeryFast"):
if (version < 7.0):
if debug: print("For QE < 7.0 RMM does not exist, reverting to davidson")
set_espresso_param(espresso_params, "diagonalization", "david")
else:
set_espresso_param(espresso_params, "diagonalization", "rmm-davidson")
elif (vasp_params[param] == "All") or (vasp_params[param] == "Conjugate"):
set_espresso_param(espresso_params, "diagonalization", "cg")
elif param == "ISPIN":
# Spin polarization
# XSpectra not compatible with nspin = 4
if (vasp_params[param] == 1):
set_espresso_param(espresso_params, "nspin", 1)
elif (vasp_params[param] == 2):
set_espresso_param(espresso_params, "nspin", 2)
elif param == "NELM":
# Max number of SCF steps
set_espresso_param(espresso_params, "electron_maxstep", vasp_params[param])
elif param == "EDIFF":
# Threshold for convergence of SCF
set_espresso_param(espresso_params, "conv_thr", vasp_params[param]/RytoeV)
elif param == "NSW":
# Number of ionic steps (calculation = scf should override this no matter what)
if vasp_params[param] == 0:
set_espresso_param(espresso_params, "nstep", 1)
else:
set_espresso_param(espresso_params, "nstep", vasp_params[param])
elif param == "ISIF":
# Calculate and print stress to outputfile
if vasp_params[param] > 0: set_espresso_param(espresso_params, "tstress", True)
elif param == "ISYM":
# Turn on/off symmetry
if vasp_params[param] < 0: set_espresso_param(espresso_params, "nosym", True)
elif param == "ENCUT":
# Cutoff for plane waves
set_espresso_param(espresso_params, "ecutwfc", round(vasp_params[param]*RytoeV, -1))
if (debug): print("Caution: cutoff on planewaves not transferable between codes, for production ecutwfc should be converged.")
elif param == "NBANDS":
# Number of bands to compute
set_espresso_param(espresso_params, "nbnd", vasp_params[param])
elif param == "ISMEAR":
# Type of smearing to apply
if (vasp_params[param] >= -1):
#espresso_params["occupations"] = "smearing"
set_espresso_param(espresso_params, "occupations", "smearing")
if vasp_params[param] > 0:
set_espresso_param(espresso_params, "smearing", "methfessel-paxton")
if (debug): print("Caution: Espresso only allows for first-order MP smearning, we have order %d" % (vasp_params[param]))
elif vasp_params[param] == 0:
set_espresso_param(espresso_params, "smearing", "gaussian")
elif vasp_params[param] == -1:
set_espresso_param(espresso_params,"smearing", "fermi-dirac")
elif (vasp_params[param] == -4):
set_espresso_param(espresso_params, "occupations", "tetrahedra_opt")
elif (vasp_params[param] == -5):
#espresso_params["occupations"] = "tetrahedra"
set_espresso_param(espresso_params,"occupations", "tetrahedra")
else:
print("Bizzare vasp smearning used, contact [email protected]")
elif param == "SIGMA":
# Value of sigma used for smearing
set_espresso_param(espresso_params, "degauss", vasp_params[param]*RytoeV)
elif param == "MAGMOM":
# Initial magnetic moments
mgmom = []
for atom in enumerate(vasp_params[param]):
mgmom.append(atom[1])
set_espresso_param(espresso_params, "starting_magnetization", mgmom)
# We ignore IBRION, ICHARG because we set calculation='scf' and use espresso default
# (atomic+random) for intial wfc and density
elif (param == "IBRION") or (param == "ICHARG"):
if debug: print("%s knowingly overidden" % (param))
elif (param == "NELMIN") or (param == "PREC") or (param == "LREAL") or (param == "KPOINT_BSE") or (param == "LPEAD") or (param == "LEFG") or (param == "LCALCPOL") or (param == "LCALCEPS") or (param == "EFIELD_PEAD"):
if debug: print("%s knowingly ignored" % (param))
elif (param == "LWAVE") or (param == "LCHARG") or (param == "LORBIT") or (param == "LAECHG"):
if debug: print("For %s, disk input/output handled by qe disk_io" % (param))
else:
print("I have ignored VASP parameter: %s" % (str(param)))
if debug: print(espresso_params)
return(espresso_params)
def set_starting_magnetization(params, atoms):
# set initial magnetization in atoms object and expunge magmom from espresso parameters dictionary
mgmom = params["starting_magnetization"]
atoms.set_initial_magnetic_moments(magmoms=mgmom)
return(atoms)
def get_vasp_kgrid(grid_data, debug=False):
# Figure out the grid definition
if grid_data["generation_style"] == "Monkhorst":
grid = grid_data["kpoints"][0]
offset = grid_data["shift"]
elif grid_data["generation_style"] == "Gamma":
grid = grid_data["kpoints"][0]
offset = np.zeros(3, dtype=int)
if debug: print(offset)
for sk in enumerate(offset):
if sk[1] == 0.5:
offset[sk[0]] = 1
if debug: print(grid, offset)
return(grid, offset)