Coverage for /builds/ase/ase/ase/calculators/dmol.py : 19.05%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""This module defines an ASE interface to DMol3.
3Contacts
4--------
5Adam Arvidsson <adam.arvidsson@chalmers.se>
6Erik Fransson <erikfr@chalmers.se>
7Anders Hellman <anders.hellman@chalmers.se>
10DMol3 environment variables
11----------------------------
12DMOL_COMMAND should point to the RunDmol script and specify the number of cores
13to prallelize over
15export DMOL_COMMAND="./RunDmol.sh -np 16"
18Example
19--------
20>>> from ase.build import bulk
21>>> from ase.calculators import DMol3
23>>> atoms = bulk('Al','fcc')
24>>> calc = DMol3()
25>>> atoms.calc = calc
26>>> print 'Potential energy %5.5f eV' % atoms.get_potential_energy()
29DMol3 calculator functionality
30-------------------------------
31This calculator does support all the functionality in DMol3.
33Firstly this calculator is limited to only handling either fully
34periodic structures (pbc = [1,1,1]) or non periodic structures (pbc=[0,0,0]).
36Internal relaxations are not supported by the calculator,
37only support for energy and forces is implemented.
39Reading eigenvalues and kpts are supported.
40Be careful with kpts and their directions (see internal coordinates below).
42Outputting the full electron density or specific bands to .grd files can be
43achieved with the plot command. The .grd files can be converted to the cube
44format using grd_to_cube().
47DMol3 internal coordinates
48---------------------------
49DMol3 may change the atomic positions / cell vectors in order to satisfy
50certain criterion ( e.g. molecule symmetry axis along z ). Specifically this
51happens when using Symmetry on/auto. This means the forces read from .grad
52will be in a different coordinates system compared to the atoms object used.
53To solve this the rotation matrix that converts the dmol coordinate system
54to the ase coordinate system is found and applied to the forces.
56For non periodic structures (pbc=[0,0,0]) the rotation matrix can be directly
57parsed from the .rot file.
58For fully periodic structures the rotation matrix is found by reading the
59cell vectors and positions used by dmol and then solving the matrix problem
60DMol_atoms * rot_mat = ase_atoms
63DMol3 files
64------------
65The supported DMol3 file formats are:
67car structure file - Angstrom and cellpar description of cell.
68incoor structure file - Bohr and cellvector describption of cell.
69 Note: incoor file not used if car file present.
70outmol outfile from DMol3 - atomic units (Bohr and Hartree)
71grad outfile for forces from DMol3 - forces in Hartree/Bohr
72grd outfile for orbitals from DMol3 - cellpar in Angstrom
74"""
76import os
77import re
78import numpy as np
79from ase import Atoms
80from ase.io import read
81from ase.io.dmol import write_dmol_car, write_dmol_incoor
82from ase.units import Hartree, Bohr
83from ase.calculators.calculator import FileIOCalculator, Parameters, ReadError
86class DMol3(FileIOCalculator):
87 """ DMol3 calculator object. """
89 implemented_properties = ['energy', 'forces']
90 default_parameters = {'functional': 'pbe',
91 'symmetry': 'on'}
92 discard_results_on_any_change = True
94 if 'DMOL_COMMAND' in os.environ:
95 command = os.environ['DMOL_COMMAND'] + ' PREFIX > PREFIX.out'
96 else:
97 command = None
99 def __init__(self, restart=None,
100 ignore_bad_restart_file=FileIOCalculator._deprecated,
101 label='dmol_calc/tmp', atoms=None, **kwargs):
102 """ Construct DMol3 calculator. """
103 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
104 label, atoms, **kwargs)
106 # tracks if DMol transformed coordinate system
107 self.internal_transformation = False
109 def write_input(self, atoms, properties=None, system_changes=None):
111 if not (np.all(atoms.pbc) or not np.any(atoms.pbc)):
112 raise RuntimeError('PBC must be all true or all false')
114 self.clean() # Remove files from old run
115 self.internal_transformation = False
116 self.ase_positions = atoms.positions.copy()
117 self.ase_cell = atoms.cell.copy()
119 FileIOCalculator.write_input(self, atoms, properties, system_changes)
121 if np.all(atoms.pbc):
122 write_dmol_incoor(self.label + '.incoor', atoms)
123 elif not np.any(atoms.pbc):
124 write_dmol_car(self.label + '.car', atoms)
126 self.write_input_file()
127 self.parameters.write(self.label + '.parameters.ase')
129 def write_input_file(self):
130 """ Writes the input file. """
131 with open(self.label + '.input', 'w') as fd:
132 self._write_input_file(fd)
134 def _write_input_file(self, fd):
135 fd.write('%-32s %s\n' % ('calculate', 'gradient'))
137 # if no key about eigs
138 fd.write('%-32s %s\n' % ('print', 'eigval_last_it'))
140 for key, value in self.parameters.items():
141 if isinstance(value, str):
142 fd.write('%-32s %s\n' % (key, value))
143 elif isinstance(value, (list, tuple)):
144 for val in value:
145 fd.write('%-32s %s\n' % (key, val))
146 else:
147 fd.write('%-32s %r\n' % (key, value))
149 def read(self, label):
150 FileIOCalculator.read(self, label)
151 geometry = self.label + '.car'
152 output = self.label + '.outmol'
153 force = self.label + '.grad'
155 for filename in [force, output, geometry]:
156 if not os.path.isfile(filename):
157 raise ReadError
159 self.atoms = read(geometry)
160 self.parameters = Parameters.read(self.label + 'parameters.ase')
161 self.read_results()
163 def read_results(self):
164 finished, message = self.finished_successfully()
165 if not finished:
166 raise RuntimeError('DMol3 run failed, see outmol file for'
167 ' more info\n\n%s' % message)
169 self.find_dmol_transformation()
170 self.read_energy()
171 self.read_forces()
173 def finished_successfully(self):
174 """ Reads outmol file and checks if job completed or failed.
176 Returns
177 -------
178 finished (bool): True if job completed, False if something went wrong
179 message (str): If job failed message contains parsed errors, else empty
181 """
182 finished = False
183 message = ""
184 for line in self._outmol_lines():
185 if line.rfind('Message: DMol3 job finished successfully') > -1:
186 finished = True
187 if line.startswith('Error'):
188 message += line
189 return finished, message
191 def find_dmol_transformation(self, tol=1e-4):
192 """Finds rotation matrix that takes us from DMol internal
193 coordinates to ase coordinates.
195 For pbc = [False, False, False] the rotation matrix is parsed from
196 the .rot file, if this file doesnt exist no rotation is needed.
198 For pbc = [True, True, True] the Dmol internal cell vectors and
199 positions are parsed and compared to self.ase_cell self.ase_positions.
200 The rotation matrix can then be found by a call to the helper
201 function find_transformation(atoms1, atoms2)
203 If a rotation matrix is needed then self.internal_transformation is
204 set to True and the rotation matrix is stored in self.rotation_matrix
206 Parameters
207 ----------
208 tol (float): tolerance for check if positions and cell are the same
209 """
211 if np.all(self.atoms.pbc): # [True, True, True]
212 dmol_atoms = self.read_atoms_from_outmol()
213 if (np.linalg.norm(self.atoms.positions - dmol_atoms.positions) <
214 tol) and (np.linalg.norm(self.atoms.cell -
215 dmol_atoms.cell) < tol):
216 self.internal_transformation = False
217 else:
218 R, err = find_transformation(dmol_atoms, self.atoms)
219 if abs(np.linalg.det(R) - 1.0) > tol:
220 raise RuntimeError('Error: transformation matrix does'
221 ' not have determinant 1.0')
222 if err < tol:
223 self.internal_transformation = True
224 self.rotation_matrix = R
225 else:
226 raise RuntimeError('Error: Could not find dmol'
227 ' coordinate transformation')
228 elif not np.any(self.atoms.pbc): # [False,False,False]
229 try:
230 data = np.loadtxt(self.label + '.rot')
231 except IOError:
232 self.internal_transformation = False
233 else:
234 self.internal_transformation = True
235 self.rotation_matrix = data[1:].transpose()
237 def read_atoms_from_outmol(self):
238 """ Reads atomic positions and cell from outmol file and returns atoms
239 object.
241 If no cell vectors are found in outmol the cell is set to np.eye(3) and
242 pbc 000.
244 Formatting for cell in outmol :
245 translation vector [a0] 1 5.1 0.0 5.1
246 translation vector [a0] 2 5.1 5.1 0.0
247 translation vector [a0] 3 0.0 5.1 5.1
249 Formatting for positions in outmol:
250 df ATOMIC COORDINATES (au)
251 df x y z
252 df Si 0.0 0.0 0.0
253 df Si 1.3 3.5 2.2
254 df binding energy -0.2309046Ha
256 Returns
257 -------
258 atoms (Atoms object): read atoms object
259 """
261 lines = self._outmol_lines()
262 found_cell = False
263 cell = np.zeros((3, 3))
264 symbols = []
265 positions = []
266 pattern_translation_vectors = re.compile(r'\s+translation\s+vector')
267 pattern_atomic_coordinates = re.compile(r'df\s+ATOMIC\s+COORDINATES')
269 for i, line in enumerate(lines):
270 if pattern_translation_vectors.match(line):
271 cell[int(line.split()[3]) - 1, :] = \
272 np.array([float(x) for x in line.split()[-3:]])
273 found_cell = True
274 if pattern_atomic_coordinates.match(line):
275 for ind, j in enumerate(range(i + 2, i + 2 + len(self.atoms))):
276 flds = lines[j].split()
277 symbols.append(flds[1])
278 positions.append(flds[2:5])
279 atoms = Atoms(symbols=symbols, positions=positions, cell=cell)
280 atoms.positions *= Bohr
281 atoms.cell *= Bohr
283 if found_cell:
284 atoms.pbc = [True, True, True]
285 atoms.wrap()
286 else:
287 atoms.pbc = [False, False, False]
288 return atoms
290 def read_energy(self):
291 """ Find and return last occurrence of Ef in outmole file. """
292 energy_regex = re.compile(r'^Ef\s+(\S+)Ha')
293 found = False
294 for line in self._outmol_lines():
295 match = energy_regex.match(line)
296 if match:
297 energy = float(match.group(1))
298 found = True
299 if not found:
300 raise RuntimeError('Could not read energy from outmol')
301 self.results['energy'] = energy * Hartree
303 def read_forces(self):
304 """ Read forces from .grad file. Applies self.rotation_matrix if
305 self.internal_transformation is True. """
306 with open(self.label + '.grad', 'r') as fd:
307 lines = fd.readlines()
309 forces = []
310 for i, line in enumerate(lines):
311 if line.startswith('$gradients'):
312 for j in range(i + 1, i + 1 + len(self.atoms)):
313 # force = - grad(Epot)
314 forces.append(np.array(
315 [-float(x) for x in lines[j].split()[1:4]]))
317 forces = np.array(forces) * Hartree / Bohr
318 if self.internal_transformation:
319 forces = np.dot(forces, self.rotation_matrix)
320 self.results['forces'] = forces
322 def get_eigenvalues(self, kpt=0, spin=0):
323 return self.read_eigenvalues(kpt, spin, 'eigenvalues')
325 def get_occupations(self, kpt=0, spin=0):
326 return self.read_eigenvalues(kpt, spin, 'occupations')
328 def get_k_point_weights(self):
329 return self.read_kpts(mode='k_point_weights')
331 def get_bz_k_points(self):
332 raise NotImplementedError
334 def get_ibz_k_points(self):
335 return self.read_kpts(mode='ibz_k_points')
337 def get_spin_polarized(self):
338 return self.read_spin_polarized()
340 def get_fermi_level(self):
341 return self.read_fermi()
343 def get_energy_contributions(self):
344 return self.read_energy_contributions()
346 def get_xc_functional(self):
347 return self.parameters['functional']
349 def read_eigenvalues(self, kpt=0, spin=0, mode='eigenvalues'):
350 """Reads eigenvalues from .outmol file.
352 This function splits into two situations:
353 1. We have no kpts just the raw eigenvalues ( Gamma point )
354 2. We have eigenvalues for each k-point
356 If calculation is spin_restricted then all eigenvalues
357 will be returned no matter what spin parameter is set to.
359 If calculation has no kpts then all eigenvalues
360 will be returned no matter what kpts parameter is set to.
362 Note DMol does usually NOT print all unoccupied eigenvalues.
363 Meaning number of eigenvalues for different kpts can vary.
364 """
366 assert mode in ['eigenvalues', 'occupations']
367 lines = self._outmol_lines()
368 pattern_kpts = re.compile(r'Eigenvalues for kvector\s+%d' % (kpt + 1))
369 for n, line in enumerate(lines):
371 # 1. We have no kpts
372 if line.split() == ['state', 'eigenvalue', 'occupation']:
373 spin_key = '+'
374 if self.get_spin_polarized():
375 if spin == 1:
376 spin_key = '-'
377 val_index = -2
378 if mode == 'occupations':
379 val_index = -1
380 values = []
381 m = n + 3
382 while True:
383 if lines[m].strip() == '':
384 break
385 flds = lines[m].split()
386 if flds[1] == spin_key:
387 values.append(float(flds[val_index]))
388 m += 1
389 return np.array(values)
391 # 2. We have kpts
392 if pattern_kpts.match(line):
393 val_index = 3
394 if self.get_spin_polarized():
395 if spin == 1:
396 val_index = 6
397 if mode == 'occupations':
398 val_index += 1
399 values = []
400 m = n + 2
401 while True:
402 if lines[m].strip() == '':
403 break
404 values.append(float(lines[m].split()[val_index]))
405 m += 1
406 return np.array(values)
407 return None
409 def _outmol_lines(self):
410 with open(self.label + '.outmol', 'r') as fd:
411 return fd.readlines()
413 def read_kpts(self, mode='ibz_k_points'):
414 """ Returns list of kpts coordinates or kpts weights. """
416 assert mode in ['ibz_k_points', 'k_point_weights']
417 lines = self._outmol_ines()
419 values = []
420 for n, line in enumerate(lines):
421 if line.startswith('Eigenvalues for kvector'):
422 if mode == 'ibz_k_points':
423 values.append([float(k_i)
424 for k_i in lines[n].split()[4:7]])
425 if mode == 'k_point_weights':
426 values.append(float(lines[n].split()[8]))
427 if values == []:
428 return None
429 return values
431 def read_spin_polarized(self):
432 """Reads, from outmol file, if calculation is spin polarized."""
434 lines = self._outmol_lines()
435 for n, line in enumerate(lines):
436 if line.rfind('Calculation is Spin_restricted') > -1:
437 return False
438 if line.rfind('Calculation is Spin_unrestricted') > -1:
439 return True
440 raise IOError('Could not read spin restriction from outmol')
442 def read_fermi(self):
443 """Reads the Fermi level.
445 Example line in outmol:
446 Fermi Energy: -0.225556 Ha -6.138 eV xyz text
447 """
448 lines = self._outmol_lines()
449 pattern_fermi = re.compile(r'Fermi Energy:\s+(\S+)\s+Ha')
450 for line in lines:
451 m = pattern_fermi.match(line)
452 if m:
453 return float(m.group(1)) * Hartree
454 return None
456 def read_energy_contributions(self):
457 """Reads the different energy contributions."""
459 lines = self._outmol_lines()
460 energies = dict()
461 for n, line in enumerate(lines):
462 if line.startswith('Energy components'):
463 m = n + 1
464 while not lines[m].strip() == '':
465 energies[lines[m].split('=')[0].strip()] = \
466 float(re.findall(
467 r"[-+]?\d*\.\d+|\d+", lines[m])[0]) * Hartree
468 m += 1
469 return energies
471 def clean(self):
472 """ Cleanup after dmol calculation
474 Only removes dmol files in self.directory,
475 does not remove the directory itself
476 """
477 file_extensions = ['basis', 'car', 'err', 'grad', 'input', 'inatm',
478 'incoor', 'kpoints', 'monitor', 'occup', 'outmol',
479 'outatom', 'rot', 'sdf', 'sym', 'tpotl', 'tpdensk',
480 'torder', 'out', 'parameters.ase']
481 files_to_clean = ['DMol3.log', 'stdouterr.txt', 'mpd.hosts']
483 files = [os.path.join(self.directory, f) for f in files_to_clean]
484 files += [''.join((self.label, '.', ext)) for ext in file_extensions]
486 for f in files:
487 try:
488 os.remove(f)
489 except OSError:
490 pass
493# Helper functions
494# ------------------
496def find_transformation(atoms1, atoms2, verbose=False, only_cell=False):
497 """ Solves Ax = B where A and B are cell and positions from atoms objects.
499 Uses numpys least square solver to solve the problem Ax = B where A and
500 B are cell vectors and positions for atoms1 and atoms2 respectively.
502 Parameters
503 ----------
504 atoms1 (Atoms object): First atoms object (A)
505 atoms2 (Atoms object): Second atoms object (B)
506 verbose (bool): If True prints for each i A[i], B[i], Ax[i]
507 only_cell (bool): If True only cell in used, otherwise cell and positions.
509 Returns
510 -------
511 x (np.array((3,3))): Least square solution to Ax = B
512 error (float): The error calculated as np.linalg.norm(Ax-b)
514 """
516 if only_cell:
517 N = 3
518 elif len(atoms1) != len(atoms2):
519 raise RuntimeError('Atoms object must be of same length')
520 else:
521 N = len(atoms1) + 3
523 # Setup matrices A and B
524 A = np.zeros((N, 3))
525 B = np.zeros((N, 3))
526 A[0:3, :] = atoms1.cell
527 B[0:3, :] = atoms2.cell
528 if not only_cell:
529 A[3:, :] = atoms1.positions
530 B[3:, :] = atoms2.positions
532 # Solve least square problem Ax = B
533 lstsq_fit = np.linalg.lstsq(A, B, rcond=-1)
534 x = lstsq_fit[0]
535 error = np.linalg.norm(np.dot(A, x) - B)
537 # Print comparison between A, B and Ax
538 if verbose:
539 print('%17s %33s %35s %24s' % ('A', 'B', 'Ax', '|Ax-b|'))
540 for a, b in zip(A, B):
541 ax = np.dot(a, x)
542 loss = np.linalg.norm(ax - b)
543 print('(', end='')
544 for a_i in a:
545 print('%8.5f' % a_i, end='')
546 print(') (', end='')
547 for b_i in b:
548 print('%8.5f ' % b_i, end='')
549 print(') (', end='')
550 for ax_i in ax:
551 print('%8.5f ' % ax_i, end='')
552 print(') %8.5f' % loss)
554 return x, error
557def grd_to_file(atoms, grd_file, new_file):
558 """ Reads grd_file and converts data to cube format and writes to
559 cube_file.
561 Note: content of grd_file and atoms object are assumed to match with the
562 same orientation.
564 Parameters
565 -----------
566 atoms (Atoms object): atoms object grd_file data is for
567 grd_file (str): filename of .grd file
568 new_file (str): filename to write grd-data to, must be ASE format
569 that supports data argument
570 """
571 from ase.io import write
573 atoms_copy = atoms.copy()
574 data, cell, origin = read_grd(grd_file)
575 atoms_copy.cell = cell
576 atoms_copy.positions += origin
577 write(new_file, atoms_copy, data=data)
580def read_grd(filename):
581 """ Reads .grd file
583 Notes
584 -----
585 origin_xyz is offset with half a grid point in all directions to be
586 compatible with the cube format
587 Periodic systems is not guaranteed to be oriented correctly
588 """
589 from ase.geometry.cell import cellpar_to_cell
591 with open(filename, 'r') as fd:
592 lines = fd.readlines()
594 cell_data = np.array([float(fld) for fld in lines[2].split()])
595 cell = cellpar_to_cell(cell_data)
596 grid = [int(fld) + 1 for fld in lines[3].split()]
597 data = np.empty(grid)
599 origin_data = [int(fld) for fld in lines[4].split()[1:]]
600 origin_xyz = cell[0] * (-float(origin_data[0]) - 0.5) / (grid[0] - 1) + \
601 cell[1] * (-float(origin_data[2]) - 0.5) / (grid[1] - 1) + \
602 cell[2] * (-float(origin_data[4]) - 0.5) / (grid[2] - 1)
604 # Fastest index describes which index ( x or y ) varies fastest
605 # 1: x , 3: y
606 fastest_index = int(lines[4].split()[0])
607 assert fastest_index in [1, 3]
608 if fastest_index == 3:
609 grid[0], grid[1] = grid[1], grid[0]
611 dummy_counter = 5
612 for i in range(grid[2]):
613 for j in range(grid[1]):
614 for k in range(grid[0]): # Fastest index
615 if fastest_index == 1:
616 data[k, j, i] = float(lines[dummy_counter])
617 elif fastest_index == 3:
618 data[j, k, i] = float(lines[dummy_counter])
619 dummy_counter += 1
621 return data, cell, origin_xyz
624if __name__ == '__main__':
625 from ase.build import molecule
627 atoms = molecule('H2')
628 calc = DMol3()
629 atoms.calc = calc
630 # ~ 60 sec calculation
631 print('Potential energy %5.5f eV' % atoms.get_potential_energy())