Coverage for /builds/debichem-team/python-ase/ase/calculators/gulp.py: 25.23%
218 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +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'
126 gin_path = os.path.join(self.directory, self.prefix + '.gin')
127 with open(gin_path, 'w') as fd:
128 fd.write(s)
130 def read_results(self):
131 FileIOCalculator.read(self, self.label)
132 got_path = os.path.join(self.directory, self.label + '.got')
133 if not os.path.isfile(got_path):
134 raise ReadError
136 with open(got_path) as fd:
137 lines = fd.readlines()
139 cycles = -1
140 self.optimized = None
141 for i, line in enumerate(lines):
142 m = re.match(r'\s*Total lattice energy\s*=\s*(\S+)\s*eV', line)
143 if m:
144 energy = float(m.group(1))
145 self.results['energy'] = energy
146 self.results['free_energy'] = energy
148 elif line.find('Optimisation achieved') != -1:
149 self.optimized = True
151 elif line.find('Final Gnorm') != -1:
152 self.Gnorm = float(line.split()[-1])
154 elif line.find('Cycle:') != -1:
155 cycles += 1
157 elif line.find('Final Cartesian derivatives') != -1:
158 s = i + 5
159 forces = []
160 while True:
161 s = s + 1
162 if lines[s].find("------------") != -1:
163 break
164 if lines[s].find(" s ") != -1:
165 continue
166 g = lines[s].split()[3:6]
167 G = [-float(x) * eV / Ang for x in g]
168 forces.append(G)
169 forces = np.array(forces)
170 self.results['forces'] = forces
172 elif line.find('Final internal derivatives') != -1:
173 s = i + 5
174 forces = []
175 while True:
176 s = s + 1
177 if lines[s].find("------------") != -1:
178 break
179 g = lines[s].split()[3:6]
181 # Uncomment the section below to separate the numbers when
182 # there is no space between them, in the case of long
183 # numbers. This prevents the code to break if numbers are
184 # too big.
186 '''for t in range(3-len(g)):
187 g.append(' ')
188 for j in range(2):
189 min_index=[i+1 for i,e in enumerate(g[j][1:])
190 if e == '-']
191 if j==0 and len(min_index) != 0:
192 if len(min_index)==1:
193 g[2]=g[1]
194 g[1]=g[0][min_index[0]:]
195 g[0]=g[0][:min_index[0]]
196 else:
197 g[2]=g[0][min_index[1]:]
198 g[1]=g[0][min_index[0]:min_index[1]]
199 g[0]=g[0][:min_index[0]]
200 break
201 if j==1 and len(min_index) != 0:
202 g[2]=g[1][min_index[0]:]
203 g[1]=g[1][:min_index[0]]'''
205 G = [-float(x) * eV / Ang for x in g]
206 forces.append(G)
207 forces = np.array(forces)
208 self.results['forces'] = forces
210 elif line.find('Final cartesian coordinates of atoms') != -1:
211 s = i + 5
212 positions = []
213 while True:
214 s = s + 1
215 if lines[s].find("------------") != -1:
216 break
217 if lines[s].find(" s ") != -1:
218 continue
219 xyz = lines[s].split()[3:6]
220 XYZ = [float(x) * Ang for x in xyz]
221 positions.append(XYZ)
222 positions = np.array(positions)
223 self.atoms.set_positions(positions)
225 elif line.find('Final stress tensor components') != -1:
226 res = [0., 0., 0., 0., 0., 0.]
227 for j in range(3):
228 var = lines[i + j + 3].split()[1]
229 res[j] = float(var)
230 var = lines[i + j + 3].split()[3]
231 res[j + 3] = float(var)
232 stress = np.array(res)
233 self.results['stress'] = stress
235 elif line.find('Final Cartesian lattice vectors') != -1:
236 lattice_vectors = np.zeros((3, 3))
237 s = i + 2
238 for j in range(s, s + 3):
239 temp = lines[j].split()
240 for k in range(3):
241 lattice_vectors[j - s][k] = float(temp[k])
242 self.atoms.set_cell(lattice_vectors)
243 if self.fractional_coordinates is not None:
244 self.fractional_coordinates = np.array(
245 self.fractional_coordinates)
246 self.atoms.set_scaled_positions(
247 self.fractional_coordinates)
249 elif line.find('Final fractional coordinates of atoms') != -1:
250 s = i + 5
251 scaled_positions = []
252 while True:
253 s = s + 1
254 if lines[s].find("------------") != -1:
255 break
256 if lines[s].find(" s ") != -1:
257 continue
258 xyz = lines[s].split()[3:6]
259 XYZ = [float(x) for x in xyz]
260 scaled_positions.append(XYZ)
261 self.fractional_coordinates = scaled_positions
263 self.steps = cycles
265 def get_opt_state(self):
266 return self.optimized
268 def get_opt_steps(self):
269 return self.steps
271 def get_Gnorm(self):
272 return self.Gnorm
274 def library_check(self):
275 if self.parameters['library'] is not None:
276 if 'GULP_LIB' not in self.cfg:
277 raise RuntimeError("Be sure to have set correctly $GULP_LIB "
278 "or to have the force field library.")
281class Conditions:
282 """Atomic labels for the GULP calculator.
284 This class manages an array similar to
285 atoms.get_chemical_symbols() via get_atoms_labels() method, but
286 with atomic labels in stead of atomic symbols. This is useful
287 when you need to use calculators like GULP or lammps that use
288 force fields. Some force fields can have different atom type for
289 the same element. In this class you can create a set_rule()
290 function that assigns labels according to structural criteria."""
292 def __init__(self, atoms):
293 self.atoms = atoms
294 self.atoms_symbols = atoms.get_chemical_symbols()
295 self.atoms_labels = atoms.get_chemical_symbols()
296 self.atom_types = []
298 def min_distance_rule(self, sym1, sym2,
299 ifcloselabel1=None, ifcloselabel2=None,
300 elselabel1=None, max_distance=3.0):
301 """Find pairs of atoms to label based on proximity.
303 This is for, e.g., the ffsioh or catlow force field, where we
304 would like to identify those O atoms that are close to H
305 atoms. For each H atoms, we must specially label one O atom.
307 This function is a rule that allows to define atom labels (like O1,
308 O2, O_H etc..) starting from element symbols of an Atoms
309 object that a force field can use and according to distance
310 parameters.
312 Example:
313 atoms = read('some_xyz_format.xyz')
314 a = Conditions(atoms)
315 a.set_min_distance_rule('O', 'H', ifcloselabel1='O2',
316 ifcloselabel2='H', elselabel1='O1')
317 new_atoms_labels = a.get_atom_labels()
319 In the example oxygens O are going to be labeled as O2 if they
320 are close to a hydrogen atom othewise are labeled O1.
322 """
324 if ifcloselabel1 is None:
325 ifcloselabel1 = sym1
326 if ifcloselabel2 is None:
327 ifcloselabel2 = sym2
328 if elselabel1 is None:
329 elselabel1 = sym1
331 # self.atom_types is a list of element types used instead of element
332 # symbols in orger to track the changes made. Take care of this because
333 # is very important.. gulp_read function that parse the output
334 # has to know which atom_type it has to associate with which
335 # atom_symbol
336 #
337 # Example: [['O','O1','O2'],['H', 'H_C', 'H_O']]
338 # this because Atoms oject accept only atoms symbols
339 self.atom_types.append([sym1, ifcloselabel1, elselabel1])
340 self.atom_types.append([sym2, ifcloselabel2])
342 dist_mat = self.atoms.get_all_distances()
343 index_assigned_sym1 = []
344 index_assigned_sym2 = []
346 for i in range(len(self.atoms_symbols)):
347 if self.atoms_symbols[i] == sym2:
348 dist_12 = 1000
349 index_assigned_sym2.append(i)
350 for t in range(len(self.atoms_symbols)):
351 if (self.atoms_symbols[t] == sym1
352 and dist_mat[i, t] < dist_12
353 and t not in index_assigned_sym1):
354 dist_12 = dist_mat[i, t]
355 closest_sym1_index = t
356 index_assigned_sym1.append(closest_sym1_index)
358 for i1, i2 in zip(index_assigned_sym1, index_assigned_sym2):
359 if dist_mat[i1, i2] > max_distance:
360 raise ValueError('Cannot unambiguously apply minimum-distance '
361 'rule because pairings are not obvious. '
362 'If you wish to ignore this, then increase '
363 'max_distance.')
365 for s in range(len(self.atoms_symbols)):
366 if s in index_assigned_sym1:
367 self.atoms_labels[s] = ifcloselabel1
368 elif s not in index_assigned_sym1 and self.atoms_symbols[s] == sym1:
369 self.atoms_labels[s] = elselabel1
370 elif s in index_assigned_sym2:
371 self.atoms_labels[s] = ifcloselabel2
373 def get_atom_types(self):
374 return self.atom_types
376 def get_atoms_labels(self):
377 labels = np.array(self.atoms_labels)
378 return labels