Commit c0990c44 authored by Ask Hjorth Larsen's avatar Ask Hjorth Larsen
Browse files

autofix flakes in calculators

parent 81b576c0
......@@ -59,7 +59,8 @@ class ACE(FileIOCalculator):
changed_parameters = FileIOCalculator.set(self, **kwargs)
# Add default values for repeated parameter sections with self.default_parameters using order.
# Also add empty dictionary as an indicator for section existence if no relevant default_parameters exist.
# Also add empty dictionary as an indicator for section existence if no
# relevant default_parameters exist.
if 'order' in kwargs:
new_parameters['order'] = kwargs['order']
section_sets = set(kwargs['order'])
......@@ -67,7 +68,8 @@ class ACE(FileIOCalculator):
repeat = kwargs['order'].count(section_name)
if section_name in self.default_parameters.keys():
for i in range(repeat-1):
new_parameters[section_name] += deepcopy(self.default_parameters[section_name])
new_parameters[section_name] += deepcopy(
self.default_parameters[section_name])
else:
new_parameters[section_name] = []
for i in range(repeat):
......@@ -81,7 +83,8 @@ class ACE(FileIOCalculator):
i = 0
for section_param in kwargs[section]:
new_parameters[section][i] = update_parameter(new_parameters[section][i], section_param)
new_parameters[section][i] = update_parameter(
new_parameters[section][i], section_param)
i += 1
self.parameters = new_parameters
return changed_parameters
......@@ -93,12 +96,14 @@ class ACE(FileIOCalculator):
with open(filename, 'r') as fd:
lines = fd.readlines()
if 'WARNING' in lines:
raise ReadError("Not convergy energy in log file {}.".format(filename))
raise ReadError(
"Not convergy energy in log file {}.".format(filename))
if '! total energy' not in lines:
raise ReadError("Wrong ACE-Molecule log file {}.".format(filename))
if not os.path.isfile(filename):
raise ReadError("Wrong ACE-Molecule input file {}.".format(filename))
raise ReadError(
"Wrong ACE-Molecule input file {}.".format(filename))
self.read_results()
......@@ -133,9 +138,11 @@ class ACE(FileIOCalculator):
Updated version of self.parameters; geometry file and optionally Force section are updated.
'''
copied_parameters = deepcopy(self.parameters)
if properties is not None and "forces" in properties and 'Force' not in copied_parameters['order']:
if properties is not None and "forces" in properties and 'Force' not in copied_parameters[
'order']:
copied_parameters['order'].append('Force')
copied_parameters["BasicInformation"][0]["GeometryFilename"] = "{}.xyz".format(self.label)
copied_parameters["BasicInformation"][0]["GeometryFilename"] = "{}.xyz".format(
self.label)
copied_parameters["BasicInformation"][0]["GeometryFormat"] = "xyz"
return copied_parameters
......@@ -150,7 +157,7 @@ class ACE(FileIOCalculator):
'''
filename = self.label + '.log'
#quantities = ['energy', 'forces', 'atoms', 'excitation-energy']
#for section_name in quantities:
# for section_name in quantities:
#self.results = read_acemolecule_out(filename)
self.results = read(filename, format='acemolecule-out')
......@@ -164,16 +171,30 @@ class ACE(FileIOCalculator):
depth: Nested input depth.
'''
for section, section_param in section.items():
if isinstance(section_param, str) or isinstance(section_param, int) or isinstance(section_param, float):
fpt.write(' ' * depth + str(section) + " " + str(section_param) + "\n")
if isinstance(section_param, str) or isinstance(
section_param, int) or isinstance(section_param, float):
fpt.write(
' ' *
depth +
str(section) +
" " +
str(section_param) +
"\n")
else:
if isinstance(section_param, dict):
fpt.write(' ' * depth + "%% " + str(section) + "\n")
self.write_acemolecule_section(fpt, section_param, depth + 1)
self.write_acemolecule_section(
fpt, section_param, depth + 1)
fpt.write(' ' * depth + "%% End\n")
if isinstance(section_param, list):
for val in section_param:
fpt.write(' ' * depth + str(section) + " " + str(val) + "\n")
fpt.write(
' ' *
depth +
str(section) +
" " +
str(val) +
"\n")
def write_acemolecule_input(self, fpt, param, depth=0):
'''Write ACE-Molecule input
......@@ -270,7 +291,8 @@ def update_parameter(oldpar, newpar):
for section, section_param in newpar.items():
if section in oldpar:
if isinstance(section_param, dict):
oldpar[section] = update_parameter(oldpar[section], section_param)
oldpar[section] = update_parameter(
oldpar[section], section_param)
else:
oldpar[section] = section_param
else:
......
......@@ -32,7 +32,7 @@ class LippincottStuttman:
'Al': 0.533,
'Si': 0.583,
}
def __call__(self, el1: str, el2: str,
length: float) -> Tuple[float, float]:
"""Bond polarizability
......@@ -106,12 +106,12 @@ class Linearized:
length0, al, ald, ap, apd = self._data[bond]
return al + ald * (length - length0), ap + apd * (length - length0)
class BondPolarizability(StaticPolarizabilityCalculator):
def __init__(self, model=LippincottStuttman()):
self.model = model
def __call__(self, atoms, radiicut=1.5):
"""Sum up the bond polarizability from all bonds
......
......@@ -621,7 +621,8 @@ class Calculator(BaseCalculator):
self.prefix = None
if label is not None:
if self.directory == '.' and '/' in label:
# We specified directory in label, and nothing in the diretory key
# We specified directory in label, and nothing in the diretory
# key
self.label = label
elif '/' not in label:
# We specified our directory in the directory keyword
......
......@@ -227,7 +227,7 @@ Keyword Description
3 = no attempt is made to look for
castep_keywords.json at all
``castep_keywords`` Can be used to pass a CastepKeywords object that is
then used with no attempt to actually load a
then used with no attempt to actually load a
castep_keywords.json file. Most useful for debugging
and testing purposes.
......@@ -782,7 +782,7 @@ End CASTEP Interface Documentation
kpts = kpts.copy()
if (kpts.get('spacing') is not None
and kpts.get('density') is not None):
and kpts.get('density') is not None):
raise ValueError(
'Cannot set kpts spacing and density simultaneously.')
else:
......@@ -1492,7 +1492,7 @@ End CASTEP Interface Documentation
# Read in eigenvalues from bands file
bands_file = castep_file[:-7] + '.bands'
if (self.param.task.value is not None
and self.param.task.value.lower() == 'bandstructure'):
and self.param.task.value.lower() == 'bandstructure'):
self._band_structure = self.band_structure(bandfile=bands_file)
else:
try:
......
......@@ -166,11 +166,11 @@ class CombineMM(Calculator):
if not self.mask[i]:
ct2 += 1
if ((ct2 == self.apm2) &
(self.apm2 != self.atoms2.calc.sites_per_mol)):
(self.apm2 != self.atoms2.calc.sites_per_mol)):
virtual_mask.append(False)
ct2 = 0
if ((ct1 == self.apm1) &
(self.apm1 != self.atoms1.calc.sites_per_mol)):
(self.apm1 != self.atoms1.calc.sites_per_mol)):
virtual_mask.append(True)
ct1 = 0
......
......@@ -541,6 +541,7 @@ class Cp2kShell:
class InputSection:
"""Represents a section of a CP2K input file"""
def __init__(self, name, params=None):
self.name = name.upper()
self.params = params
......
......@@ -189,7 +189,7 @@ class PureDFTD3(FileIOCalculator):
# stresses) can be bypassed. This will greatly speed up calculations
# in dense 3D-periodic systems with three-body corrections. But, we
# can no longer say that we implement forces and stresses.
#if not self.parameters['grad']:
# if not self.parameters['grad']:
# for val in ['forces', 'stress']:
# if val in self.implemented_properties:
# self.implemented_properties.remove(val)
......
......@@ -948,7 +948,8 @@ End EAM Interface Documentation
plt.plot(curvex, curvey[i](curvex), label=label)
plt.legend()
def multielem_subplot(self, curvex, curvey, xlabel, ylabel, name, plt, half=True):
def multielem_subplot(self, curvex, curvey, xlabel,
ylabel, name, plt, half=True):
plt.xlabel(xlabel)
plt.ylabel(ylabel)
for i in np.arange(self.Nelements):
......
......@@ -5,6 +5,7 @@ from ase.units import Hartree, Bohr
class Excitation:
"""Base class for a single excitation"""
def __init__(self, energy, index, mur, muv=None, magn=None):
"""
Parameters
......@@ -65,7 +66,7 @@ class Excitation:
magn = np.array([float(l.pop(0)) for i in range(3)])
except IndexError:
magn = None
return cls(energy, index, mur, muv, magn)
def get_dipole_me(self, form='r'):
......@@ -92,10 +93,11 @@ class Excitation:
class ExcitationList(list):
"""Base class for excitions from the ground state"""
def __init__(self):
# initialise empty list
super().__init__()
# set default energy scale to get eV
self.energy_to_eV_scale = 1.
......
......@@ -53,7 +53,7 @@ class Exciting:
def update(self, atoms):
if (not self.converged or
len(self.numbers) != len(atoms) or
(self.numbers != atoms.get_atomic_numbers()).any()):
(self.numbers != atoms.get_atomic_numbers()).any()):
self.initialize(atoms)
self.calculate(atoms)
elif ((self.positions != atoms.get_positions()).any() or
......@@ -96,9 +96,9 @@ class Exciting:
assert (Path(self.dir) / 'INFO.OUT').is_file()
assert (Path(self.dir) / 'info.xml').exists()
#syscall = ('cd %(dir)s; %(bin)s;' %
# syscall = ('cd %(dir)s; %(bin)s;' %
# {'dir': self.dir, 'bin': self.excitingbinary})
#print(syscall)
# print(syscall)
#assert os.system(syscall) == 0
self.read()
......
......@@ -49,6 +49,7 @@ class FLEUR:
(atoms to optimize etc.) are specified by hand in *inp*
"""
def __init__(self, xc='LDA', kpts=None, nbands=None, convergence=None,
width=None, kmax=None, mixer=None, maxiter=None,
maxrelax=20, workdir=None, equivatoms=True, rmt=None,
......@@ -150,22 +151,29 @@ class FLEUR:
out = p.stdout.read()
err = p.stderr.read()
print(mode, ': stat= ', stat, ' out= ', out, ' err=', err)
# special handling of exit status from density generation and regular fleur.x
# special handling of exit status from density generation and regular
# fleur.x
if mode in ['density']:
if '!' in err:
os.chdir(self.start_dir)
raise RuntimeError(executable_use + ' exited with a code %s' % err)
raise RuntimeError(
executable_use +
' exited with a code %s' %
err)
else:
if stat != 0:
os.chdir(self.start_dir)
raise RuntimeError(executable_use + ' exited with a code %d' % stat)
raise RuntimeError(
executable_use +
' exited with a code %d' %
stat)
def update(self, atoms):
"""Update a FLEUR calculation."""
if (not self.converged or
len(self.numbers) != len(atoms) or
(self.numbers != atoms.get_atomic_numbers()).any()):
(self.numbers != atoms.get_atomic_numbers()).any()):
self.initialize(atoms)
self.calculate(atoms)
elif ((self.positions != atoms.get_positions()).any() or
......@@ -277,7 +285,9 @@ class FLEUR:
while not self.converged:
if self.niter > self.itmax:
os.chdir(self.start_dir)
raise RuntimeError('FLEUR failed to convergence in %d iterations' % self.itmax)
raise RuntimeError(
'FLEUR failed to convergence in %d iterations' %
self.itmax)
self.run_executable(mode='fleur', executable='FLEUR')
......@@ -322,7 +332,9 @@ class FLEUR:
os.system('cp cdn1 cdn1_%d' % nrelax)
if nrelax > self.maxrelax:
os.chdir(self.start_dir)
raise RuntimeError('Failed to relax in %d iterations' % self.maxrelax)
raise RuntimeError(
'Failed to relax in %d iterations' %
self.maxrelax)
self.converged = False
def write_inp(self, atoms):
......@@ -376,8 +388,10 @@ class FLEUR:
else:
# generate inequivalent atoms, by using non-integer Z
# (only the integer part will be used as Z of the atom)
# see http://www.flapw.de/pm/index.php?n=User-Documentation.InputFileForTheInputGenerator
fh.write('%3d.%04d' % (Z, n)) # MDTMP don't think one can calculate more that 10**4 atoms
# see
# http://www.flapw.de/pm/index.php?n=User-Documentation.InputFileForTheInputGenerator
# MDTMP don't think one can calculate more that 10**4 atoms
fh.write('%3d.%04d' % (Z, n))
for el in pos:
fh.write(' %21.16f' % el)
fh.write('\n')
......@@ -411,7 +425,9 @@ class FLEUR:
lines[ln] = 'mjw non-relativic\n'
del lines[ln+1]
else:
raise RuntimeError('XC-functional %s is not supported' % self.xc)
raise RuntimeError(
'XC-functional %s is not supported' %
self.xc)
if line.startswith('Window'):
# few things are set around this line
window_ln = ln
......@@ -428,7 +444,8 @@ class FLEUR:
# gmaxxc cutoff for PW-expansion of XC-potential ( > 2*kmax, < gmax)
if self.kmax and line.startswith('vchk'):
gmax = 3. * self.kmax
line = ' %10.6f %10.6f\n' % (gmax, int(2.5 * self.kmax * 10)/10.)
line = ' %10.6f %10.6f\n' % (
gmax, int(2.5 * self.kmax * 10)/10.)
lines[ln-1] = line
# Fermi width
if self.width and line.startswith('gauss'):
......@@ -441,8 +458,10 @@ class FLEUR:
self.kpts[2])
lines[ln] = line
# itmax
if self.itmax < self.itmax_step_default and line.startswith('itmax'):
# decrease number of SCF steps; increasing is done by 'while not self.converged:'
if self.itmax < self.itmax_step_default and line.startswith(
'itmax'):
# decrease number of SCF steps; increasing is done by 'while not
# self.converged:'
lsplit = line.split(',')
if lsplit[0].find('itmax') != -1:
lsplit[0] = 'itmax=' + ('%2d' % self.itmax)
......
......@@ -3,6 +3,7 @@ from ase.calculators.calculator import PropertyNotImplementedError
class Calculator:
"Deprecated!!!!"
def __init__(self):
return
......
......@@ -279,7 +279,7 @@ class Gromacs(FileIOCalculator):
generate topology (self.label+'top')
and structure file in .g96 format (self.label + '.g96')
"""
#generate structure and topology files
# generate structure and topology files
# In case of predefinded topology file this is not done
subcmd = 'pdb2gmx'
command = ' '.join([
......@@ -302,7 +302,7 @@ class Gromacs(FileIOCalculator):
resulting file is self.label + '.tpr
"""
#generate gromacs run input file (gromacs.tpr)
# generate gromacs run input file (gromacs.tpr)
try:
os.remove(self.label + '.tpr')
except OSError:
......@@ -349,7 +349,7 @@ class Gromacs(FileIOCalculator):
"""Write input parameters to input file."""
FileIOCalculator.write_input(self, atoms, properties, system_changes)
#print self.parameters
# print self.parameters
with open(self.label + '.mdp', 'w') as myfile:
for key, val in self.parameters.items():
if val is not None:
......@@ -361,7 +361,7 @@ class Gromacs(FileIOCalculator):
""" set atoms and do the calculation """
# performs an update of the atoms
self.atoms = atoms.copy()
#must be g96 format for accuracy, alternatively binary formats
# must be g96 format for accuracy, alternatively binary formats
write_gromos(self.label + '.g96', atoms)
# does run to get forces and energies
self.calculate()
......@@ -405,7 +405,7 @@ class Gromacs(FileIOCalculator):
with open(self.label + '.Energy.xvg') as fd:
lastline = fd.readlines()[-1]
energy = float(lastline.split()[1])
#We go for ASE units !
# We go for ASE units !
#self.energy = energy * units.kJ / units.mol
self.results['energy'] = energy * units.kJ / units.mol
# energies are about 100 times bigger in Gromacs units
......@@ -423,7 +423,7 @@ class Gromacs(FileIOCalculator):
with open(self.label + '.Force.xvg', 'r') as fd:
lastline = fd.readlines()[-1]
forces = np.array([float(f) for f in lastline.split()[1:]])
#We go for ASE units !gromacsForce.xvg
# We go for ASE units !gromacsForce.xvg
#self.forces = np.array(forces)/ units.nm * units.kJ / units.mol
#self.forces = np.reshape(self.forces, (-1, 3))
tmp_forces = forces / units.nm * units.kJ / units.mol
......
......@@ -168,7 +168,10 @@ class GULP(FileIOCalculator):
break
g = lines[s].split()[3:6]
# Uncomment the section below to separate the numbers when there is no space between them, in the case of long numbers. This prevents the code to break if numbers are too big.
# Uncomment the section below to separate the numbers when
# there is no space between them, in the case of long
# numbers. This prevents the code to break if numbers are
# too big.
'''for t in range(3-len(g)):
g.append(' ')
......@@ -227,7 +230,8 @@ class GULP(FileIOCalculator):
lattice_vectors[j-s][k] = float(temp[k])
self.atoms.set_cell(lattice_vectors)
if self.fractional_coordinates is not None:
self.fractional_coordinates = np.array(self.fractional_coordinates)
self.fractional_coordinates = np.array(
self.fractional_coordinates)
self.atoms.set_scaled_positions(self.fractional_coordinates)
elif line.find('Final fractional coordinates of atoms') != -1:
......@@ -312,8 +316,8 @@ class Conditions:
if elselabel1 is None:
elselabel1 = sym1
#self.atom_types is a list of element types used instead of element
#symbols in orger to track the changes made. Take care of this because
# self.atom_types is a list of element types used instead of element
# symbols in orger to track the changes made. Take care of this because
# is very important.. gulp_read function that parse the output
# has to know which atom_type it has to associate with which
# atom_symbol
......@@ -334,7 +338,7 @@ class Conditions:
for t in range(len(self.atoms_symbols)):
if (self.atoms_symbols[t] == sym1
and dist_mat[i, t] < dist_12
and t not in index_assigned_sym1):
and t not in index_assigned_sym1):
dist_12 = dist_mat[i, t]
closest_sym1_index = t
index_assigned_sym1.append(closest_sym1_index)
......
......@@ -108,6 +108,7 @@ class H2MorseCalculator(MorsePotential):
class H2MorseExcitedStatesCalculator():
"""First singlet excited states of H2 from Morse potentials"""
def __init__(self, nstates=3):
"""
Parameters
......@@ -151,6 +152,7 @@ class H2MorseExcitedStatesCalculator():
class H2MorseExcitedStates(ExcitationList):
"""First singlet excited states of H2"""
def __init__(self, nstates=3):
"""
Parameters
......@@ -198,6 +200,7 @@ class H2Excitation(Excitation):
class H2MorseExcitedStatesAndCalculator(
H2MorseExcitedStatesCalculator, H2MorseExcitedStates):
"""Traditional joined object for backward compatibility only"""
def __init__(self, calculator, nstates=3):
if isinstance(calculator, str):
exlist = H2MorseExcitedStates.read(calculator, nstates)
......
......@@ -59,7 +59,7 @@ class HarmonicForceField:
Together with the :class:`HarmonicCalculator` this
:class:`HarmonicForceField` can be used to compute
Anharmonic Corrections to the Harmonic Approximation. [1]_
Parameters
----------
ref_atoms: :class:`~ase.Atoms` object
......@@ -182,7 +182,8 @@ class HarmonicForceField:
and back. Relevant literature:
* Peng, C. et al. J. Comput. Chem. 1996, 17 (1), 49-56.
* Baker, J. et al. J. Chem. Phys. 1996, 105 (1), 192–212."""
jac0 = self.get_jacobian(self.parameters['ref_atoms']) # Jacobian (dq/dx)
jac0 = self.get_jacobian(
self.parameters['ref_atoms']) # Jacobian (dq/dx)
jac0 = self.constrain_jac(jac0) # for reference Cartesian coordinates
ijac0 = self.get_ijac(jac0, self.parameters['rcond'])
self.transform2reference_hessians(jac0, ijac0) # perform projection
......@@ -232,12 +233,12 @@ class HarmonicForceField:
"""Return a tuple with energy and forces in Cartesian coordinates for
a given :class:`~ase.Atoms` object."""
q = self.get_q_from_x(atoms).ravel()
if self.parameters['cartesian']:
x = atoms.get_positions().ravel()
x0 = self.x0
hessian_x = self._hessian_x
if self.parameters['variable_orientation']:
# determine x0 for present orientation
x0 = self.back_transform(x, q, self.q0, atoms.copy())
......@@ -249,11 +250,11 @@ class HarmonicForceField:
self.check_redundancy(jac0) # check for coordinate failure
# determine hessian_x for present orientation
hessian_x = jac0.T @ self._hessian_q @ jac0
xdiff = x - x0
forces_x = -hessian_x @ xdiff
energy = -0.5 * (forces_x * xdiff).sum()
else:
jac = self.get_jacobian(atoms)
self.check_redundancy(jac) # check for coordinate failure
......@@ -261,7 +262,7 @@ class HarmonicForceField:
forces_q = -self._hessian_q @ qdiff
forces_x = forces_q @ jac
energy = -0.5 * (forces_q * qdiff).sum()
energy += self.parameters['ref_energy']
forces_x = forces_x.reshape(int(forces_x.size / 3), 3)
return energy, forces_x
......@@ -363,7 +364,9 @@ class SpringCalculator(Calculator):
F = 0.0
masses, counts = np.unique(self.atoms.get_masses(), return_counts=True)
for m, c in zip(masses, counts):
F += c * SpringCalculator.compute_Einstein_solid_free_energy(self.k, m, T, method)
F += c * \
SpringCalculator.compute_Einstein_solid_free_energy(
self.k, m, T, method)
return F
@staticmethod
......@@ -398,7 +401,8 @@ class SpringCalculator(Calculator):
omega = np.sqrt(k / m) # angular frequency 1/s
if method == 'classical':
F_einstein = 3 * units.kB * T * np.log(hbar * omega / (units.kB * T))
F_einstein = 3 * units.kB * T * \
np.log(hbar * omega / (units.kB * T))
elif method == 'QM':
log_factor = np.log(1.0 - np.exp(-hbar * omega / (units.kB * T)))
F_einstein = 3 * units.kB * T * log_factor + 1.5 * hbar * omega
......
......@@ -218,7 +218,8 @@ class LennardJones(Calculator):
pairwise_forces = -24 * epsilon * (2 * c12 - c6) / r2 # du_ij
if smooth:
# order matters, otherwise the pairwise energy is already modified
# order matters, otherwise the pairwise energy is already
# modified
pairwise_forces = (
cutoff_fn * pairwise_forces + 2 * d_cutoff_fn * pairwise_energies
)
......@@ -238,7 +239,8 @@ class LennardJones(Calculator):
# no lattice, no stress
if self.atoms.cell.rank == 3:
stresses = full_3x3_to_voigt_6_stress(stresses)
self.results['stress'] = stresses.sum(axis=0) / self.atoms.get_volume()
self.results['stress'] = stresses.sum(
axis=0) / self.atoms.get_volume()
self.results['stresses'] = stresses / self.atoms.get_volume()
energy = energies.sum()
......@@ -268,7 +270,8 @@ def cutoff_function(r, rc, ro):
return np.where(
r < ro,
1.0,
np.where(r < rc, (rc - r) ** 2 * (rc + 2 * r - 3 * ro) / (rc - ro) ** 3, 0.0),
np.where(r < rc, (rc - r) ** 2 * (rc + 2 *