/
structure_handling.py
258 lines (224 loc) · 9.68 KB
/
structure_handling.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
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
###############################
# This file is part of PyLaDa.
#
# Copyright (C) 2013 National Renewable Energy Lab
#
# PyLaDa is a high throughput computational platform for Physics. It aims to make it easier to
# submit large numbers of jobs on supercomputers. It provides a python interface to physical input,
# such as crystal structures, as well as to a number of DFT (VASP, CRYSTAL) and atomic potential
# programs. It is able to organise and launch computational jobs on PBS and SLURM.
#
# PyLaDa is free software: you can redistribute it and/or modify it under the terms of the GNU
# General Public License as published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# PyLaDa is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
# the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along with PyLaDa. If not, see
# <http://www.gnu.org/licenses/>.
###############################
# -*- coding: utf-8 -*-
""" Reads and writes structure to and from espresso input """
__docformat__ = "restructuredtext en"
__all__ = ['read_structure', 'add_structure']
from pylada import logger
logger = logger.getChild('espresso')
""" Logger for espresso """
def read_structure(filename):
""" Reads crystal structure from Espresso input """
from numpy import dot, array
from quantities import bohr_radius as a0, angstrom, Ry
from ..crystal import Structure
from .. import error
from . import Namelist
from .card import read_cards
namelists = Namelist()
namelists.read(filename)
cards = read_cards(filename)
if 'atomic_positions' not in set([u.name for u in cards]):
raise error.RuntimeError("Card ATOMIC_POSITIONS is missing from input")
positions = [u for u in cards if u.name == 'atomic_positions'][0]
if 'cell_parameters' not in set([u.name for u in cards]):
cell_parameters = None
if namelists.system.ibrav == 0:
raise error.RuntimeError("Card CELL_PARAMETERS is missing")
else:
cell_parameters = [u for u in cards if u.name == 'cell_parameters'][0]
cell, scale = read_cell_and_scale(namelists.system, cell_parameters)
result = Structure()
result.cell = cell
result.scale = scale
for line in positions.value.rstrip().lstrip().split('\n'):
line = line.split()
result.add_atom(array(line[1:4], dtype='float64'), line[0])
if positions.subtitle == 'bohr':
factor = float(a0.rescale(result.scale))
for atom in result:
atom.pos *= factor
elif positions.subtitle == 'angstrom':
factor = float(angstrom.rescale(result.scale))
for atom in result:
atom.pos *= factor
elif positions.subtitle == 'crystal':
for atom in result:
atom.pos = dot(result.cell, atom.pos)
elif positions.subtitle == 'crystal_sg':
raise error.RuntimeError("Reading symmetric atomic positions is not implemented")
if 'atomic_forces' in set([u.name for u in cards]):
atomic_forces = [u for u in cards if u.name == 'atomic_forces'][0]
if atomic_forces.value is None:
raise error.RuntimeError("Atomic forces card is empty")
lines = atomic_forces.value.rstrip().lstrip().split('\n')
if len(lines) != len(result):
raise error.RuntimeError("Number forces and number of atoms do not match")
for atom, force in zip(result, lines):
atom.force = array(force.rstrip().lstrip().split()[1:4], dtype='float64') * Ry / a0
return result
def read_cell_and_scale(system, cell_parameters):
from numpy import identity, array
if system.ibrav == 0:
return _read_free(system, cell_parameters)
elif system.ibrav == 1:
return identity(3, dtype='float64'), _get_scale(system)
elif system.ibrav == 2:
return 0.5 * array([[-1, 0, 1], [0, 1, 1], [-1, 1, 0]], dtype='float64').transpose(),\
_get_scale(system)
elif system.ibrav == 3:
return 0.5 * array([[1, 1, 1], [-1, 1, 1], [-1, -1, 1]], dtype='float64').transpose(),\
_get_scale(system)
elif system.ibrav == 4:
return _read_hexa(system)
elif system.ibrav == 5:
logger.warning('ibrav=5 has not been tested')
return _read_trig(system)
elif system.ibrav == -5:
logger.warning('ibrav=-5 has not been tested')
return _read_mtrig(system)
else:
NotImplementedError("Reading from this kind of lattice has not been implemented")
def _get_scale(system):
from ..misc import Sequence
from quantities import bohr_radius, angstrom
celldm = getattr(system, 'celldm', None)
if celldm == 'A':
return angstrom
elif celldm is not None:
if not isinstance(celldm, Sequence):
return celldm * bohr_radius
elif len(celldm) > 0:
return celldm[0] * bohr_radius
elif hasattr(system, 'A'):
return system.a * angstrom
return bohr_radius
def _get_params(system):
""" Returns celldim(2-6), whatever the input may be """
from numpy import zeros
from quantities import angstrom
result = zeros(5, dtype='float64')
celldim = getattr(system, 'celldm', None)
if celldim is not None:
result[:len(celldim) - 1] = celldim[1:]
return result
scale = _get_scale(system)
result[0] = float((getattr(system, 'b', 0) * angstrom / scale).simplified)
result[1] = float((getattr(system, 'c', 0) * angstrom / scale).simplified)
result[2] = getattr(system, 'cosab', 0)
result[3] = getattr(system, 'cosac', 0)
result[4] = getattr(system, 'cosbc', 0)
return result
def _read_free(system, cell_parameters):
""" Read free cell """
from numpy import array
from quantities import bohr_radius, angstrom
scale = _get_scale(system)
if cell_parameters.subtitle is not None:
if cell_parameters.subtitle == 'bohr':
scale = bohr_radius
elif cell_parameters.subtitle == 'angstrom':
scale = angstrom
cell = array([u.split() for u in cell_parameters.value.rstrip().lstrip().split('\n')],
dtype='float64').transpose()
return cell, scale
def _read_hexa(system):
from numpy import sqrt, array
scale = _get_scale(system)
c_over_a = _get_params(system)[1]
cell = array([[1, 0, 0], [-0.5, sqrt(3.) / 2., 0], [0, 0, c_over_a]], dtype='float64')
return cell.transpose(), scale
def _read_trig(system):
from numpy import sqrt, array
scale = _get_scale(system)
c = _get_params(system)[2]
tx, ty, tz = sqrt((1. - c) / 2.), sqrt((1. - c) / 6.), sqrt((1 + 2 * c) / 3.)
cell = array([[tx, -ty, tz], [0, 2. * ty, tz], [-tx, -ty, tz]], dtype='float64')
return cell.transpose(), scale
def _read_mtrig(system):
from numpy import sqrt, array
scale = _get_scale(system)
c = _get_params(system)[2]
ty, tz = sqrt((1. - c) / 6.), sqrt((1 + 2 * c) / 3.)
u, v = tz - 2 * sqrt(2) * ty, tz + sqrt(2) * ty
cell = array([[u, v, v], [v, u, v], [v, v, u]], dtype='float64')
return cell.transpose(), scale / sqrt(3.)
def add_structure(structure, f90namelist, cards):
""" Modifies f90namelist and cards according to structure """
from . import Card
from f90nml import Namelist as F90Namelist
from quantities import bohr_radius
if 'system' not in f90namelist:
f90namelist['system'] = F90Namelist()
for key in ['a', 'b', 'c', 'cosab', 'cosac', 'cosbc', 'celldm', 'nat', 'ntyp']:
f90namelist['system'].pop(key, None)
f90namelist['system'].pop(key.upper(), None)
f90namelist['system']['ibrav'] = 0
f90namelist['system']['celldm'] = float(structure.scale.rescale(bohr_radius))
f90namelist['system']['nat'] = len(structure)
f90namelist['system']['ntyp'] = len(set([u.type for u in structure]))
card_dict = {card.name: card for card in cards}
if 'cell' not in card_dict:
cell = Card('cell_parameters')
cards.append(cell)
else:
cell = card_dict['cell']
cell.subtitle = 'alat'
cell.value = ""
for i in range(structure.cell.shape[1]):
cell.value += " ".join([str(u) for u in structure.cell[:, i]]) + "\n"
positions = card_dict.get('atomic_positions', Card('atomic_positions'))
if 'atomic_positions' not in card_dict:
positions = Card('atomic_positions')
cards.append(positions)
else:
positions = card_dict['atomic_positions']
positions.subtitle = 'alat'
positions.value = ""
for atom in structure:
positions.value += "%s %18.12e %18.12e %18.12e\n" % (atom.type, atom.pos[0], atom.pos[1],
atom.pos[2])
__add_forces_to_input(cards, structure)
def __add_forces_to_input(cards, structure):
from numpy import allclose
from quantities import Ry, bohr_radius as a0
from .card import Card
if structure is None:
return
forces = []
units = Ry / a0
for atom in structure:
force = getattr(atom, 'force', [0, 0, 0])
if hasattr(force, 'rescale'):
force = force.rescale(units).magnitude
forces.append(force)
if allclose(forces, 0):
return
atomic_forces = Card('atomic_forces', value="")
for atom, force in zip(structure, forces):
atomic_forces.value += "%s %18.15e %18.15e %18.15e\n" % (
atom.type, force[0], force[1], force[2])
# filter cards in-place: we need to modify the input sequence itself
for i, u in enumerate(list(cards)):
if u.name in 'atomic_forces':
cards.pop(i)
cards.append(atomic_forces)