Commit 991adea3 authored by Jens Jørgen Mortensen's avatar Jens Jørgen Mortensen

Merge branch 'E501' into 'master'

E501

See merge request !601
parents bfe4fe54 5ceaca88
Pipeline #98986600 failed with stage
in 3 minutes and 36 seconds
......@@ -7,10 +7,12 @@ master:
- pip install git+https://gitlab.com/ase/[email protected] > ase-install.output
- mkdir ~/.gpaw
- echo "scalapack = True; libraries += ['scalapack-openmpi']" > ~/.gpaw/siteconfig.py
- python setup.py install # > gpaw-install.output
- python setup.py install > gpaw-install.output
- useradd -m user
- su user -c 'gpaw install-data --register gpaw-datasets'
- su user -c 'gpaw -P 1 info'
- su user -c 'gpaw test --range linalg/gemm_complex.py,lcao/dos.py'
- flake8 --ignore E722,E303,E221,E731,E203,E262,E202,E402,E502,E127,E201,E305,E241,W291,E302,E225,E128,E261,E251,E501,E265,E226,E231,E129,E741,W293,W503,W504 --exclude "doc/platforms/*" gpaw doc
- flake8 --exit-zero --exclude "doc/platforms/*,gpaw/lrtddft2/*" gpaw doc > f8.out
- cat f8.out | wc -l | tee nerr
- python -c "assert int(open('nerr').read()) <= 2361, 'Please run flake8 on your code'"
- python -We:invalid -m compileall -f -q gpaw/
......@@ -17,7 +17,7 @@ class Box:
self.owns = []
self.position = None
def set_position(self, position):
self.position = np.asarray(position)
......@@ -46,7 +46,7 @@ class MPL:
def plot(self):
a4 = 100 * np.array([2**-1.75, 2**-2.25])
inch = 2.54
self.fig = plt.figure(1, a4 / inch)
self.ax = ax = self.fig.add_axes([0, 0, 1, 1], frameon=False)
ax.set_xlim(0, a4[0])
......@@ -84,7 +84,7 @@ class MPL:
bbox = b.text.get_window_extent()
t = b.text.get_transform()
b.size = t.inverted().transform(bbox.size)
for b in self.boxes:
for other, name, s, style in b.owns:
d = other.position - b.position
......@@ -116,6 +116,7 @@ def box(*args, **kwargs):
boxes.append(b)
return b
atoms = box('Atoms', [''], ['positions, numbers, cell, pbc'],
color='white')
paw = box('PAW', [], [], 'green')
......@@ -153,7 +154,10 @@ fd = box('FDWaveFunctions',
[], 'magenta')
pt = box('LFC', r'$\tilde{p}_i^a(\mathbf{r})$', [], 'red')
lcao = box('LCAOWaveFunctions',
r"$\tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r})=\sum_{\mu\mathbf{R}} C_{\sigma\mathbf{k}n\mu} \Phi_\mu(\mathbf{r} - \mathbf{R}) \exp(i\mathbf{k}\cdot\mathbf{R})$",
r"$\tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r})=" +
r"\sum_{\mu\mathbf{R}} " +
r"C_{\sigma\mathbf{k}n\mu} \Phi_\mu(\mathbf{r} - \mathbf{R}) " +
r"\exp(i\mathbf{k}\cdot\mathbf{R})$",
['S_qMM, T_qMM, P_aqMi'], 'magenta')
atoms0 = box('Atoms', '(copy)', ['positions, numbers, cell, pbc'],
color='grey')
......@@ -162,8 +166,9 @@ forces = box('ForceCalculator')
occupations = box(
'OccupationNumbers',
r'$\epsilon_{\sigma\mathbf{k}n} \rightarrow f_{\sigma\mathbf{k}n}$')
poisson = box('PoissonSolver',
r'$\nabla^2 \tilde{v}_H(\mathbf{r}) = -4\pi \tilde{\rho}(\mathbf{r})$')
poisson = box(
'PoissonSolver',
r'$\nabla^2 \tilde{v}_H(\mathbf{r}) = -4\pi \tilde{\rho}(\mathbf{r})$')
eigensolver = box('EigenSolver')
symmetry = box('Symmetry')
restrictor = box('Transformer', '(fine -> coarse)',
......@@ -175,7 +180,7 @@ kin = box('FDOperator', r'$-\frac{1}{2}\nabla^2$')
hsoperator = box('HSOperator',
[r"$\langle \psi_n | A | \psi_{n'} \rangle$",
r"$\sum_{n'}U_{nn'}|\tilde{\psi}_{n'}\rangle$"])
overlap = box('Overlap')
basisfunctions = box('BasisFunctions', r'$\Phi_\mu(\mathbf{r})$',
color='red')
......
......@@ -6,6 +6,7 @@ machine = os.environ.get('MACHINE', 'TEST')
ncores = os.environ.get('NCORES', 8)
if rank == 0:
os.chdir(machine)
os.system('python ../memory_bandwidth.py --runs=5 --startcores='+str(ncores))
#os.system('python ../memory_bandwidth.py --runs=5') # full memory benchmark
shutil.copy('memory_bandwidth_'+machine+'_py.png', '..')
os.system(
'python ../memory_bandwidth.py --runs=5 --startcores=' +
str(ncores))
shutil.copy('memory_bandwidth_' + machine + '_py.png', '..')
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import print_function
import Numeric as num
from gpaw import GPAW
from gpaw.spherical_harmonics import Y
......@@ -60,10 +59,12 @@ for n in paw.nuclei:
p.plot(X + x, num.dot(P_i, phi_i), C + '-', lw=2,
label=r'$\psi^%s$' % s.symbol)
C = 'r'
p.plot([-d / 2], [0], 'go', ms=2*0.8*4*80/a*1.0*0.53, mfc=None, label='_nolegend_')
p.plot([d / 2], [0], 'ro', ms=2*0.8*4*80/a*1.2*0.53, mfc=None, label='_nolegend_')
p.plot([-d / 2], [0], 'go', ms=2*0.8*4*80/a*1.0*0.53, mfc=None,
label='_nolegend_')
p.plot([d / 2], [0], 'ro', ms=2*0.8*4*80/a*1.2*0.53, mfc=None,
label='_nolegend_')
p.legend(loc='best')
p.xlabel(u'x [Å]')
p.ylabel(r'$\psi$')
#p.show()
# p.show()
p.savefig('co_wavefunctions.png', dpi=dpi)
# Creates: cl_field.ind_Ffe.png, qm_field.ind_Ffe.png, tot_field.ind_Ffe.png
# -*- coding: utf-8 -*-
from gpaw.mpi import world
assert world.size == 1, 'This script should be run in serial mode (with one process).'
import numpy as np
import matplotlib.pyplot as plt
......@@ -9,6 +8,8 @@ import matplotlib.pyplot as plt
from gpaw.inducedfield.inducedfield_base import BaseInducedField
from gpaw.tddft.units import aufrequency_to_eV
assert world.size == 1, 'This script should be run in serial mode.'
# Helper function
def do_plot(d_g, ng, box, atoms):
......@@ -38,6 +39,7 @@ def do_plot(d_g, ng, box, atoms):
plt.ylim([y[0], y[-1]])
ax.set_aspect('equal')
for fname, name in zip(['cl_field.ind', 'qm_field.ind', 'tot_field.ind'],
['Classical subsystem', 'Quantum subsystem',
'Total hybrid system']):
......
# Creates: field.ind_Ffe.png
# -*- coding: utf-8 -*-
from gpaw.mpi import world
assert world.size == 1, 'This script should be run in serial mode (with one process).'
import numpy as np
import matplotlib.pyplot as plt
......@@ -9,6 +8,8 @@ import matplotlib.pyplot as plt
from gpaw.inducedfield.inducedfield_base import BaseInducedField
from gpaw.tddft.units import aufrequency_to_eV
assert world.size == 1, 'This script should be run in serial mode.'
# Helper function
def do_plot(d_g, ng, box, atoms):
......@@ -38,6 +39,7 @@ def do_plot(d_g, ng, box, atoms):
plt.ylim([y[0], y[-1]])
ax.set_aspect('equal')
for fname, name in zip(['field.ind'], ['Classical system']):
# Read InducedField object
ind = BaseInducedField(fname, readmode='all')
......
......@@ -48,4 +48,6 @@ f2 = lr.lr_transitions.get_transition_contributions(0)
for (ip, val) in enumerate(f2):
if (val > 1e-3):
parprint(' %5d => %5d %lf %%\n ' %
(lr.ks_singles.kss_list[ip].occ_ind, lr.ks_singles.kss_list[ip].unocc_ind, val / 2. * 100))
(lr.ks_singles.kss_list[ip].occ_ind,
lr.ks_singles.kss_list[ip].unocc_ind,
val / 2. * 100))
......@@ -7,7 +7,9 @@ manually instead of using gpawtransport, which currently does not work
from ase import Atoms
from gpaw import GPAW, Mixer, FermiDirac
from gpaw.lcao.tools import remove_pbc, get_lcao_hamiltonian, get_lead_lcao_hamiltonian
from gpaw.lcao.tools import (remove_pbc,
get_lcao_hamiltonian,
get_lead_lcao_hamiltonian)
import pickle as pickle
a = 2.41 # Pt binding length
......@@ -48,7 +50,7 @@ H -= Ef * S
remove_pbc(atoms, H, S, 0)
# Dump the Hamiltonian and Scattering matrix to a pickle file
pickle.dump((H.astype(complex), S.astype(complex)),
pickle.dump((H.astype(complex), S.astype(complex)),
open('scat_hs.pickle', 'wb'), 2)
########################
......
......@@ -6,7 +6,8 @@ Lattice="4.12605070262 0.0 0.0 -2.06302535131 3.57326472577 0.0 0.0 0.0 18.13127
V 0.00000000 -0.00000000 9.06563772 1.00000000 23 0.00000000 0.00000000 0.00000003 2.63358805 -0.00000000
I 2.06302535 1.19108824 10.63303664 0.10000000 53 0.00000000 0.00000000 0.00792029 -0.03990743 0.00000000
I -0.00000000 2.38217648 7.49823879 0.10000000 53 -0.00000000 0.00000000 -0.00792046 -0.03990743 -0.00000000
""")
""") # noqa
Path('CrI3.xyz').write_text("""\
8
Lattice="7.00794410138 5.69308307302e-19 0.0 -3.50397205069 6.0690576201 0.0 -7.55612427101e-19 0.0 18.0635365764" Properties=species:S:1:pos:R:3:initial_magmoms:R:1:Z:I:1:forces:R:3:magmoms:R:1 energy=-31.3985425600531 dipole="-0.00239032749259 0.00138005621976 3.42821467339e-07" free_energy=-31.3985446012773 pbc="T T F" magmom=6.00000000876689 stress="-6.01056204961e-05 0.0 -6.48091502334e-21 0.0 -6.01056204961e-05 -2.59236600933e-20 -6.48091502334e-21 -2.59236600933e-20 6.05488486747e-05"
......@@ -18,4 +19,4 @@ I -2.24281648 3.88469895 7.46488009 1.00000000 53
I 2.52228800 -0.00001162 10.59865633 1.00000000 53 -0.00469075 0.00014771 0.00214167 -0.05825936
I 2.24281799 3.88469794 10.59865633 1.00000000 53 0.00247329 0.00398846 0.00214167 -0.05825936
I -1.26113394 2.18437130 10.59865633 1.00000000 53 0.00221746 -0.00413616 0.00214167 -0.05825936
""")
""") # noqa
......@@ -7,12 +7,13 @@ from gpaw.solvation import (
SolvationGPAW, # the solvation calculator
EffectivePotentialCavity, # cavity using an effective potential
Power12Potential, # a specific effective potential
LinearDielectric, # rule to construct permittivity function from the cavity
GradientSurface, # rule to calculate the surface area from the cavity
SurfaceInteraction # rule to calculate non-electrostatic interactions
LinearDielectric, # rule to construct permittivity func from the cavity
GradientSurface, # rule to calculate the surface area from the cavity
SurfaceInteraction # rule to calculate non-electrostatic interactions
)
# all parameters on the user side of the solvation API follow the ASE unit conventions (eV, Angstrom, ...)
# all parameters on the user side of the solvation API follow the ASE
# unit conventions (eV, Angstrom, ...)
# non-solvent related DFT parameters
h = 0.2
......@@ -21,14 +22,20 @@ vac = 5.0
# solvent parameters for water from J. Chem. Phys. 141, 174108 (2014)
u0 = 0.180 # eV
epsinf = 78.36 # dimensionless
gamma = 18.4 * 1e-3 * Pascal * m # convert from dyne / cm to eV / Angstrom ** 2
gamma = 18.4 * 1e-3 * Pascal * m # convert from dyne / cm to eV / Angstrom**2
T = 298.15 # Kelvin
vdw_radii = vdw_radii.copy()
vdw_radii[1] = 1.09
# atomic_radii expected by gpaw.solvation.Power12Potential have to be a callable
# mapping an Atoms object to an iterable of floats representing the atomic radii of
# the atoms in the same order as in the Atoms object (in Angstrom)
atomic_radii = lambda atoms: [vdw_radii[n] for n in atoms.numbers]
# atomic_radii expected by gpaw.solvation.Power12Potential have to be a
# callable mapping an Atoms object to an iterable of floats representing
# the atomic radii of the atoms in the same order as in the Atoms object
# (in Angstrom)
def atomic_radii(atoms):
return [vdw_radii[n] for n in atoms.numbers]
# create Atoms object for ethanol and add vacuum
atoms = molecule('CH3CH2OH')
......@@ -38,7 +45,8 @@ atoms.center(vacuum=vac)
atoms.calc = GPAW(xc='PBE', h=h, txt='gasphase.txt')
Egasphase = atoms.get_potential_energy()
# perform calculation with continuum solvent model from J. Chem. Phys. 141, 174108 (2014)
# perform calculation with continuum solvent model from
# J. Chem. Phys. 141, 174108 (2014)
atoms.calc = SolvationGPAW(
xc='PBE', h=h, txt='water.txt',
cavity=EffectivePotentialCavity(
......
......@@ -3,14 +3,17 @@ from gpaw import GPAW
from gpaw.response.df import DielectricFunction
# Part 1: Ground state calculation
atoms = bulk('Si', 'diamond', a=5.431) # Generate diamond crystal structure for silicon
calc = GPAW(mode='pw', kpts=(4,4,4)) # GPAW calculator initialization
atoms = bulk('Si', 'diamond', a=5.431)
calc = GPAW(mode='pw', kpts=(4, 4, 4))
atoms.set_calculator(calc)
atoms.get_potential_energy() # Ground state calculation is performed
calc.write('si.gpw', 'all') # Use 'all' option to write wavefunction
atoms.get_potential_energy() # Ground state calculation is performed
calc.write('si.gpw', 'all') # Use 'all' option to write wavefunction
# Part 2 : Spectrum calculation # DF: dielectric function object
df = DielectricFunction(calc='si.gpw', # Ground state gpw file (with wavefunction) as input
# Part 2 : Spectrum calculation
# DF: dielectric function object
# Ground state gpw file (with wavefunction) as input
df = DielectricFunction(calc='si.gpw',
domega0=0.05) # Using nonlinear frequency grid
df.get_dielectric_function() # By default, a file called 'df.csv' is generated
# By default, a file called 'df.csv' is generated
df.get_dielectric_function()
......@@ -3,14 +3,16 @@ import numpy as np
import pickle
import matplotlib.pyplot as plt
result = pickle.load(open('BN_GW0_results.pckl','rb'))
result = pickle.load(open('BN_GW0_results.pckl', 'rb'))
nite = result['iqp'].shape[0]
gap = []
for i in range(nite):
gap = np.append(gap,np.min(result['iqp'][i,0,:,3])-np.max(result['iqp'][i,0,:,2]))
gap = np.append(gap,
np.min(result['iqp'][i, 0, :, 3]) -
np.max(result['iqp'][i, 0, :, 2]))
plt.figure()
plt.plot(range(nite),gap,'ko-')
plt.plot(range(nite), gap, 'ko-')
plt.ylabel('Band gap (eV)')
plt.xlabel('Iteration')
plt.savefig('BN_GW0.png')
......
# Creates: gaps.csv
from __future__ import print_function
from ase.io import read
......@@ -9,6 +8,7 @@ for name in ['no_u', 'normalized_u', 'not_normalized_u']:
n = read(name + '.txt')
gap = n.calc.get_eigenvalues(spin=1)[1] - n.calc.get_eigenvalues(spin=0)[1]
gaps.append(gap)
print('%s, %.3f' % (name.replace('_', ' ').replace('u', 'U'), gap), file=fd)
print('%s, %.3f' % (name.replace('_', ' ').replace('u', 'U'), gap),
file=fd)
assert abs(gaps[1] - gaps[0] - 6.0) < 0.8
......@@ -5,7 +5,8 @@ from gpaw.xc.rpa import RPACorrelation
f = paropen('con_freq.dat', 'w')
for N in [4, 6, 8, 12, 16, 24, 32]:
rpa = RPACorrelation('N2.gpw', txt='rpa_N2_frequencies.txt', nfrequencies=N)
rpa = RPACorrelation('N2.gpw', txt='rpa_N2_frequencies.txt',
nfrequencies=N)
E = rpa.calculate(ecut=[50])
print(N, E[0], file=f)
if N == 16:
......
from __future__ import print_function
import numpy as np
from ase.units import Bohr, Hartree
from ase.parallel import paropen
......@@ -79,7 +77,8 @@ class ExteriorElectronDensity:
'relative EED weight: eed_weight']), end=' ', file=out)
print(
'#; n k s weight energy occ eed_weight', file=out)
'#; n k s weight energy occ eed_weight',
file=out)
for kpt in wfs.kpt_u:
for n in range(wfs.bd.nbands):
print('%4d %3d %1d %8.5f %10.5f %10.5f %10.5f' %
......
from __future__ import print_function
from math import sqrt
import numpy as np
......@@ -196,7 +195,7 @@ class SimpleStm(STM):
"""
return self.scan_const_density(self.current_to_density(current),
bias)
def scan_const_density(self, density, bias):
"""Get the height image for constant density [e/Angstrom^3].
"""
......@@ -204,7 +203,7 @@ class SimpleStm(STM):
self.heights = heights / Bohr
return self.heights
def write(self, file=None):
"""Write STM data to a file in gnuplot readable tyle."""
......@@ -237,7 +236,8 @@ class SimpleStm(STM):
print('#', file=f)
print('# density=', self.density, '[e/Angstrom^3]', end=' ', file=f)
print(
'(current=', self.density_to_current(self.density), '[nA])', file=f)
'(current=', self.density_to_current(self.density), '[nA])',
file=f)
print('# x[Angs.] y[Angs.] h[Angs.] (-1 is not found)', file=f)
for i in range(nx):
for j in range(ny):
......
......@@ -72,7 +72,8 @@ def vdWradii(symbols, xc):
n += 1
# linear interpolation
ncut = (n_g[n - 1] +
(n_g[n] - n_g[n - 1]) * (R - r_g[n - 1]) / (r_g[n] - r_g[n - 1]))
(n_g[n] - n_g[n - 1]) * (R - r_g[n - 1]) /
(r_g[n] - r_g[n - 1]))
# print "Z, Zrg, ncut", Z, Zrg, ncut
# find own R at this density
......@@ -82,7 +83,8 @@ def vdWradii(symbols, xc):
n += 1
# linear interpolation
R = (r_g[n - 1] +
(r_g[n] - r_g[n - 1]) * (ncut - n_g[n - 1]) / (n_g[n] - n_g[n - 1]))
(r_g[n] - r_g[n - 1]) * (ncut - n_g[n - 1]) /
(n_g[n] - n_g[n - 1]))
radius[symbol] = R * Bohr
radii.append(radius[symbol])
......
......@@ -77,15 +77,9 @@ class AllElectron:
f = 1.0
else:
f = float(conf[2:])
#try:
assert n == self.n_j[j]
assert l == self.l_j[j]
self.f_j[j] = f
#except IndexError:
# self.n_j.append(n)
# self.l_j.append(l)
# self.f_j.append(f)
# self.e_j.append(self.e_j[-1])
j += 1
else:
j += {'He': 1,
......@@ -107,7 +101,6 @@ class AllElectron:
self.f_j = [self.Z]
self.e_j = [self.e_j[-1]]
t = self.text
t()
if scalarrel:
......@@ -251,7 +244,7 @@ class AllElectron:
mix = 0.01
nitermax = 2000
e_j[0] /= self.tw_coeff
if Z > 10 : #help convergence for third row elements
if Z > 10: # help convergence for third row elements
mix = 0.002
nitermax = 10000
......@@ -342,11 +335,12 @@ class AllElectron:
if self.orbital_free:
# e and vr are not scaled back
# instead Ekin is scaled for total energy (printed and inside setup)
# instead Ekin is scaled for total energy
# (printed and inside setup)
Ekin *= self.tw_coeff
t()
t('Lambda:{0}'.format(self.tw_coeff))
t('Correct eigenvalue:{0}'.format(e_j[0]*self.tw_coeff))
t('Correct eigenvalue:{0}'.format(e_j[0] * self.tw_coeff))
t()
t()
......@@ -891,6 +885,7 @@ guess for the density).
return nodes, A
if __name__ == '__main__':
a = AllElectron('Cu', scalarrel=True)
a.run()
......@@ -103,7 +103,8 @@ def get_orbitals_by_energy_shift(opts, setup, **kwargs):
phit_g[1:] /= rgd.r_g[1:]
gcut = rgd.ceil(rsplit)
#tailnorm = np.dot(rgd.dr_g[gcut:],
# (rgd.r_g[gcut:] * orbital_bf.phit_g[gcut:])**2)**0.5
# (rgd.r_g[gcut:] *
# orbital_bf.phit_g[gcut:])**2)**0.5
#print 'tailnorm', tailnorm
dphit_g = orbital_bf.phit_g[:gcut+1] - phit_g[:gcut+1]
bf = BasisFunction(l=orbital_bf.l,
......
# flake8: noqa
import functools
from ase.calculators.calculator import Calculator
......
# flake8: noqa
'''A class for computing the
parameters for Marcus theory
from two constrained DFT
......
......@@ -15,7 +15,6 @@ from ase.io import Trajectory
from gpaw import GPAW
from gpaw.mpi import serial_comm, rank
from gpaw.kpt_descriptor import KPointDescriptor
#from gpaw.dfpt.poisson import PoissonSolver, FFTPoissonSolver
from gpaw.poisson import PoissonSolver, FFTPoissonSolver
from gpaw.dfpt.responsecalculator import ResponseCalculator
from gpaw.dfpt.phononperturbation import PhononPerturbation
......@@ -25,6 +24,7 @@ from gpaw.dfpt.electronphononcoupling import ElectronPhononCoupling
from gpaw.symmetry import Symmetry
class PhononCalculator:
"""This class defines the interface for phonon calculations."""
......@@ -294,8 +294,10 @@ class PhononCalculator:
def get_filename_string(self):
"""Return string template for force constant filenames."""
name_str = (self.name + '.' + 'q_%%0%ii_' % len(str(self.kd.nibzkpts)) +
'a_%%0%ii_' % len(str(len(self.atoms))) + 'v_%i' + '.pckl')
name_str = (self.name + '.' +
'q_%%0%ii_' % len(str(self.kd.nibzkpts)) +
'a_%%0%ii_' % len(str(len(self.atoms))) +
'v_%i' + '.pckl')
return name_str
......@@ -416,7 +418,7 @@ class PhononCalculator:
# multiply with mass prefactor
u_nav = u_avn[:, omega2_n.argsort()].T.copy() * m_inv_av
# Multiply with mass prefactor
u_kn.append(u_nav.reshape((3*N, -1, 3)))
u_kn.append(u_nav.reshape((3 * N, -1, 3)))
else:
omega2_n = la.eigvalsh(D_q, UPLO='L')
......@@ -448,7 +450,7 @@ class PhononCalculator:
return omega_kn
def write_modes(self, q_c, branches=0, kT=units.kB*300, repeat=(1, 1, 1),
def write_modes(self, q_c, branches=0, kT=units.kB * 300, repeat=(1, 1, 1),
nimages=30, acoustic=True):
"""Write mode to trajectory file.
......@@ -511,11 +513,12 @@ class PhononCalculator:
mode_av = np.zeros((len(self.atoms), 3), dtype=self.dtype)
indices = self.dyn.get_indices()
mode_av[indices] = u_av
mode_mav = (np.vstack([mode_av]*M) * phase_ma[:, np.newaxis]).real
mode_mav = (np.vstack([mode_av] * M) *
phase_ma[:, np.newaxis]).real
traj = Trajectory('%s.mode.%d.traj' % (self.name, n), 'w')
for x in np.linspace(0, 2*pi, nimages, endpoint=False):
for x in np.linspace(0, 2 * pi, nimages, endpoint=False):
# XXX Is it correct to take out the sine component here ?
atoms.set_positions(pos_mav + sin(x) * mode_mav)
traj.write(atoms)
......
......@@ -24,8 +24,8 @@ class PhononPerturbation(Perturbation):
def __init__(self, calc, kd, poisson_solver, dtype=float, **kwargs):
"""Store useful objects, e.g. lfc's for the various atomic functions.
Depending on whether the system is periodic or finite, Poisson's equation
is solved with FFT or multigrid techniques, respectively.
Depending on whether the system is periodic or finite, Poisson's
equation is solved with FFT or multigrid techniques, respectively.
Parameters
----------
......
......@@ -66,7 +66,8 @@ class Phonons(phonons.Phonons):
"""Read force constants from files."""
# Data structure for force constants
self.C_qaavv = [dict([(a, dict([(a_, np.zeros((3, 3), dtype=self.dtype))
self.C_qaavv = [dict([(a, dict([(a_, np.zeros((3, 3),
dtype=self.dtype))
for a_ in self.indices]))
for a in self.indices])
for q in range(self.kd.nibzkpts)]
......@@ -187,7 +188,8 @@ class Phonons(phonons.Phonons):
# Reshape before Fourier transforming
shape = self.D_k.shape
Dq_lmn = self.D_k.reshape(N_c + shape[1:])
DR_lmn = fft.ifftn(fft.ifftshift(Dq_lmn, axes=(0, 1, 2)), axes=(0, 1, 2))
DR_lmn = fft.ifftn(fft.ifftshift(Dq_lmn, axes=(0, 1, 2)),
axes=(0, 1, 2))
if debug:
# Check that D_R is real enough
......
from __future__ import print_function
"""This module implements a linear response calculator class."""
__all__ = ["ResponseCalculator"]
......@@ -380,15 +379,17 @@ class ResponseCalculator:
rhs_G = -1 * rhs_nG[n]
# Solve Sternheimer equation
iter, info = self.linear_solver.solve(self.sternheimer_operator,
psit1_G, rhs_G)
iter, info = self.linear_solver.solve(
self.sternheimer_operator,
psit1_G, rhs_G)
if verbose:
print("\tBand %2.1i -" % n, end=' ')
if info == 0:
if verbose:
print("linear solver converged in %i iterations" % iter)
print("linear solver converged in %i iterations" %
iter)
elif info > 0:
assert False, ("linear solver did not converge in maximum "
"number (=%i) of iterations for "
......
......@@ -118,7 +118,8 @@ class ELF:
# For periodic boundary conditions
if self.paw.wfs.kd.symmetry is not None:
self.paw.wfs.kd.symmetry.symmetrize(self.taut_sG[0], self.paw.wfs.gd)
self.paw.wfs.kd.symmetry.symmetrize(self.taut_sG[0],
self.paw.wfs.gd)
self.nt_grad2_sG[:] = 0.0
......
......@@ -200,6 +200,7 @@ class PointChargePotential(ExternalPotential):
return F_pv * Hartree / Bohr
class CDFTPotential(ExternalPotential):
# Dummy class to make cDFT compatible with new external potential class ClassName(object):
# Dummy class to make cDFT compatible with new external
# potential class ClassName(object):
def __init__(self):
self.name = 'CDFTPotential'
......@@ -250,7 +250,8 @@ class FDTDPoissonSolver:
def estimate_memory(self, mem):
# self.cl.poisson_solver.estimate_memory(mem)