Coverage for /builds/kinetik161/ase/ase/calculators/gulp.py: 25.00%
216 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""This module defines an ASE interface to GULP.
3Written by:
5Andy Cuko <andi.cuko@upmc.fr>
6Antoni Macia <tonimacia@gmail.com>
8EXPORT ASE_GULP_COMMAND="/path/to/gulp < PREFIX.gin > PREFIX.got"
10Keywords
11Options
13"""
14import os
15import re
17import numpy as np
19from ase.calculators.calculator import FileIOCalculator, ReadError
20from ase.units import Ang, eV
23class GULPOptimizer:
24 def __init__(self, atoms, calc):
25 self.atoms = atoms
26 self.calc = calc
28 def todict(self):
29 return {'type': 'optimization',
30 'optimizer': 'GULPOptimizer'}
32 def run(self, fmax=None, steps=None, **gulp_kwargs):
33 if fmax is not None:
34 gulp_kwargs['gmax'] = fmax
35 if steps is not None:
36 gulp_kwargs['maxcyc'] = steps
38 self.calc.set(**gulp_kwargs)
39 self.atoms.calc = self.calc
40 self.atoms.get_potential_energy()
41 self.atoms.cell = self.calc.get_atoms().cell
42 self.atoms.positions[:] = self.calc.get_atoms().positions
45class GULP(FileIOCalculator):
46 implemented_properties = ['energy', 'free_energy', 'forces', 'stress']
47 _legacy_default_command = 'gulp < PREFIX.gin > PREFIX.got'
48 discard_results_on_any_change = True
49 default_parameters = dict(
50 keywords='conp gradients',
51 options=[],
52 shel=[],
53 library="ffsioh.lib",
54 conditions=None)
56 def get_optimizer(self, atoms):
57 gulp_keywords = self.parameters.keywords.split()
58 if 'opti' not in gulp_keywords:
59 raise ValueError('Can only create optimizer from GULP calculator '
60 'with "opti" keyword. Current keywords: {}'
61 .format(gulp_keywords))
63 opt = GULPOptimizer(atoms, self)
64 return opt
66 def __init__(self, restart=None,
67 ignore_bad_restart_file=FileIOCalculator._deprecated,
68 label='gulp', atoms=None, optimized=None,
69 Gnorm=1000.0, steps=1000, conditions=None, **kwargs):
70 """Construct GULP-calculator object."""
71 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
72 label, atoms, **kwargs)
73 self.optimized = optimized
74 self.Gnorm = Gnorm
75 self.steps = steps
76 self.conditions = conditions
77 self.library_check()
78 self.atom_types = []
80 # GULP prints the fractional coordinates before the Final
81 # lattice vectors so they need to be stored and then atoms
82 # positions need to be set after we get the Final lattice
83 # vectors
84 self.fractional_coordinates = None
86 def write_input(self, atoms, properties=None, system_changes=None):
87 FileIOCalculator.write_input(self, atoms, properties, system_changes)
88 p = self.parameters
90 # Build string to hold .gin input file:
91 s = p.keywords
92 s += '\ntitle\nASE calculation\nend\n\n'
94 if all(self.atoms.pbc):
95 cell_params = self.atoms.cell.cellpar()
96 # Formating is necessary since Gulp max-line-length restriction
97 s += 'cell\n{:9.6f} {:9.6f} {:9.6f} ' \
98 '{:8.5f} {:8.5f} {:8.5f}\n'.format(*cell_params)
99 s += 'frac\n'
100 coords = self.atoms.get_scaled_positions()
101 else:
102 s += 'cart\n'
103 coords = self.atoms.get_positions()
105 if self.conditions is not None:
106 c = self.conditions
107 labels = c.get_atoms_labels()
108 self.atom_types = c.get_atom_types()
109 else:
110 labels = self.atoms.get_chemical_symbols()
112 for xyz, symbol in zip(coords, labels):
113 s += ' {:2} core' \
114 ' {:10.7f} {:10.7f} {:10.7f}\n' .format(symbol, *xyz)
115 if symbol in p.shel:
116 s += ' {:2} shel' \
117 ' {:10.7f} {:10.7f} {:10.7f}\n' .format(symbol, *xyz)
119 if p.library:
120 s += f'\nlibrary {p.library}\n'
122 if p.options:
123 for t in p.options:
124 s += f'{t}\n'
125 with open(self.prefix + '.gin', 'w') as fd:
126 fd.write(s)
128 def read_results(self):
129 FileIOCalculator.read(self, self.label)
130 if not os.path.isfile(self.label + '.got'):
131 raise ReadError
133 with open(self.label + '.got') as fd:
134 lines = fd.readlines()
136 cycles = -1
137 self.optimized = None
138 for i, line in enumerate(lines):
139 m = re.match(r'\s*Total lattice energy\s*=\s*(\S+)\s*eV', line)
140 if m:
141 energy = float(m.group(1))
142 self.results['energy'] = energy
143 self.results['free_energy'] = energy
145 elif line.find('Optimisation achieved') != -1:
146 self.optimized = True
148 elif line.find('Final Gnorm') != -1:
149 self.Gnorm = float(line.split()[-1])
151 elif line.find('Cycle:') != -1:
152 cycles += 1
154 elif line.find('Final Cartesian derivatives') != -1:
155 s = i + 5
156 forces = []
157 while True:
158 s = s + 1
159 if lines[s].find("------------") != -1:
160 break
161 if lines[s].find(" s ") != -1:
162 continue
163 g = lines[s].split()[3:6]
164 G = [-float(x) * eV / Ang for x in g]
165 forces.append(G)
166 forces = np.array(forces)
167 self.results['forces'] = forces
169 elif line.find('Final internal derivatives') != -1:
170 s = i + 5
171 forces = []
172 while True:
173 s = s + 1
174 if lines[s].find("------------") != -1:
175 break
176 g = lines[s].split()[3:6]
178 # Uncomment the section below to separate the numbers when
179 # there is no space between them, in the case of long
180 # numbers. This prevents the code to break if numbers are
181 # too big.
183 '''for t in range(3-len(g)):
184 g.append(' ')
185 for j in range(2):
186 min_index=[i+1 for i,e in enumerate(g[j][1:])
187 if e == '-']
188 if j==0 and len(min_index) != 0:
189 if len(min_index)==1:
190 g[2]=g[1]
191 g[1]=g[0][min_index[0]:]
192 g[0]=g[0][:min_index[0]]
193 else:
194 g[2]=g[0][min_index[1]:]
195 g[1]=g[0][min_index[0]:min_index[1]]
196 g[0]=g[0][:min_index[0]]
197 break
198 if j==1 and len(min_index) != 0:
199 g[2]=g[1][min_index[0]:]
200 g[1]=g[1][:min_index[0]]'''
202 G = [-float(x) * eV / Ang for x in g]
203 forces.append(G)
204 forces = np.array(forces)
205 self.results['forces'] = forces
207 elif line.find('Final cartesian coordinates of atoms') != -1:
208 s = i + 5
209 positions = []
210 while True:
211 s = s + 1
212 if lines[s].find("------------") != -1:
213 break
214 if lines[s].find(" s ") != -1:
215 continue
216 xyz = lines[s].split()[3:6]
217 XYZ = [float(x) * Ang for x in xyz]
218 positions.append(XYZ)
219 positions = np.array(positions)
220 self.atoms.set_positions(positions)
222 elif line.find('Final stress tensor components') != -1:
223 res = [0., 0., 0., 0., 0., 0.]
224 for j in range(3):
225 var = lines[i + j + 3].split()[1]
226 res[j] = float(var)
227 var = lines[i + j + 3].split()[3]
228 res[j + 3] = float(var)
229 stress = np.array(res)
230 self.results['stress'] = stress
232 elif line.find('Final Cartesian lattice vectors') != -1:
233 lattice_vectors = np.zeros((3, 3))
234 s = i + 2
235 for j in range(s, s + 3):
236 temp = lines[j].split()
237 for k in range(3):
238 lattice_vectors[j - s][k] = float(temp[k])
239 self.atoms.set_cell(lattice_vectors)
240 if self.fractional_coordinates is not None:
241 self.fractional_coordinates = np.array(
242 self.fractional_coordinates)
243 self.atoms.set_scaled_positions(
244 self.fractional_coordinates)
246 elif line.find('Final fractional coordinates of atoms') != -1:
247 s = i + 5
248 scaled_positions = []
249 while True:
250 s = s + 1
251 if lines[s].find("------------") != -1:
252 break
253 if lines[s].find(" s ") != -1:
254 continue
255 xyz = lines[s].split()[3:6]
256 XYZ = [float(x) for x in xyz]
257 scaled_positions.append(XYZ)
258 self.fractional_coordinates = scaled_positions
260 self.steps = cycles
262 def get_opt_state(self):
263 return self.optimized
265 def get_opt_steps(self):
266 return self.steps
268 def get_Gnorm(self):
269 return self.Gnorm
271 def library_check(self):
272 if self.parameters['library'] is not None:
273 if 'GULP_LIB' not in self.cfg:
274 raise RuntimeError("Be sure to have set correctly $GULP_LIB "
275 "or to have the force field library.")
278class Conditions:
279 """Atomic labels for the GULP calculator.
281 This class manages an array similar to
282 atoms.get_chemical_symbols() via get_atoms_labels() method, but
283 with atomic labels in stead of atomic symbols. This is useful
284 when you need to use calculators like GULP or lammps that use
285 force fields. Some force fields can have different atom type for
286 the same element. In this class you can create a set_rule()
287 function that assigns labels according to structural criteria."""
289 def __init__(self, atoms):
290 self.atoms = atoms
291 self.atoms_symbols = atoms.get_chemical_symbols()
292 self.atoms_labels = atoms.get_chemical_symbols()
293 self.atom_types = []
295 def min_distance_rule(self, sym1, sym2,
296 ifcloselabel1=None, ifcloselabel2=None,
297 elselabel1=None, max_distance=3.0):
298 """Find pairs of atoms to label based on proximity.
300 This is for, e.g., the ffsioh or catlow force field, where we
301 would like to identify those O atoms that are close to H
302 atoms. For each H atoms, we must specially label one O atom.
304 This function is a rule that allows to define atom labels (like O1,
305 O2, O_H etc..) starting from element symbols of an Atoms
306 object that a force field can use and according to distance
307 parameters.
309 Example:
310 atoms = read('some_xyz_format.xyz')
311 a = Conditions(atoms)
312 a.set_min_distance_rule('O', 'H', ifcloselabel1='O2',
313 ifcloselabel2='H', elselabel1='O1')
314 new_atoms_labels = a.get_atom_labels()
316 In the example oxygens O are going to be labeled as O2 if they
317 are close to a hydrogen atom othewise are labeled O1.
319 """
321 if ifcloselabel1 is None:
322 ifcloselabel1 = sym1
323 if ifcloselabel2 is None:
324 ifcloselabel2 = sym2
325 if elselabel1 is None:
326 elselabel1 = sym1
328 # self.atom_types is a list of element types used instead of element
329 # symbols in orger to track the changes made. Take care of this because
330 # is very important.. gulp_read function that parse the output
331 # has to know which atom_type it has to associate with which
332 # atom_symbol
333 #
334 # Example: [['O','O1','O2'],['H', 'H_C', 'H_O']]
335 # this because Atoms oject accept only atoms symbols
336 self.atom_types.append([sym1, ifcloselabel1, elselabel1])
337 self.atom_types.append([sym2, ifcloselabel2])
339 dist_mat = self.atoms.get_all_distances()
340 index_assigned_sym1 = []
341 index_assigned_sym2 = []
343 for i in range(len(self.atoms_symbols)):
344 if self.atoms_symbols[i] == sym2:
345 dist_12 = 1000
346 index_assigned_sym2.append(i)
347 for t in range(len(self.atoms_symbols)):
348 if (self.atoms_symbols[t] == sym1
349 and dist_mat[i, t] < dist_12
350 and t not in index_assigned_sym1):
351 dist_12 = dist_mat[i, t]
352 closest_sym1_index = t
353 index_assigned_sym1.append(closest_sym1_index)
355 for i1, i2 in zip(index_assigned_sym1, index_assigned_sym2):
356 if dist_mat[i1, i2] > max_distance:
357 raise ValueError('Cannot unambiguously apply minimum-distance '
358 'rule because pairings are not obvious. '
359 'If you wish to ignore this, then increase '
360 'max_distance.')
362 for s in range(len(self.atoms_symbols)):
363 if s in index_assigned_sym1:
364 self.atoms_labels[s] = ifcloselabel1
365 elif s not in index_assigned_sym1 and self.atoms_symbols[s] == sym1:
366 self.atoms_labels[s] = elselabel1
367 elif s in index_assigned_sym2:
368 self.atoms_labels[s] = ifcloselabel2
370 def get_atom_types(self):
371 return self.atom_types
373 def get_atoms_labels(self):
374 labels = np.array(self.atoms_labels)
375 return labels