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