Commit 4ee8210c authored by marc barbarosa's avatar marc barbarosa
Browse files

Merge branch 'master' of https://gitlab.com/ase_mbpt_lcao/ase

parents 97d16c2b e5457b6a
Pipeline #18763852 failed with stage
in 46 minutes and 56 seconds
Fork of ASE containing tools to plot and analyse data from the mbpt_lcao (http://mbpt-domiprod.wikidot.com/) program.
Atomic Simulation Environment
=============================
......
......@@ -1930,7 +1930,9 @@ def symbols2numbers(symbols):
symbols = string2symbols(symbols)
numbers = []
for s in symbols:
if isinstance(s, basestring):
if s =='GH':
numbers.append(0)
elif isinstance(s, basestring):
numbers.append(atomic_numbers[s])
else:
numbers.append(s)
......
......@@ -402,3 +402,4 @@ def readWFSX(fname):
return WFSX_tuple(nkpoints, nspin, norbitals, gamma, orb2atm,
orb2strspecies, orb2ao, orb2n, orb2strsym,
kpoints, DFT_E, DFT_X, mo_spin_kpoint_2_is_read)
from __future__ import division
import numpy as np
import os
from ase.utils import basestring
class MBPT_LCAO:
......@@ -43,7 +42,7 @@ class MBPT_LCAO:
for k, v in self.param.items():
if isinstance(v, np.ndarray):
f.write(k + ' {0} {1} {2}\n'.format(v[0], v[1], v[2]))
elif isinstance(v, basestring):
elif isinstance(v, str):
f.write(k + ' ' + v + '\n')
elif k == 'group_species' or k == 'species_iter':
gp = '{'
......
This diff is collapsed.
......@@ -3,6 +3,18 @@ import numpy as np
import re
def get_lattice_param(atm, guess):
lc = []
for i, pos in enumerate(atm.get_positions()):
for j, pos2 in enumerate(atm.get_positions()):
r = np.sqrt(np.dot(pos-pos2, pos-pos2))
if r<guess:
lc.append(r)
lc = np.array(lc)
return np.sum(lc)/lc.shape[0]
def read_file(fname):
"""
read the file fname and return a list of the lines.
......
from __future__ import division
import numpy as np
from ase.mbpt_lcao_utils.utils import pts_coord_to_Array_ind, Array_ind_to_pts_coord
from ase.mbpt_lcao_utils.Read_data import read_data
#
#
#
def trajectory_elec(args, data, dt):
"""
Calculate the electron trajectory incide the box
"""
t = 0
R0 = t*args.tem_input['vnorm']*args.tem_input['vdir'] + args.tem_input['b']
while (R0[0]>data.box[0, 0] and R0[1]>data.box[0, 1] and R0[2]>data.box[0, 2]):
t = t - dt
R0 = t*args.tem_input['vnorm']*args.tem_input['vdir'] + args.tem_input['b']
tmin = t+dt
t = 0
R0 = t*args.tem_input['vnorm']*args.tem_input['vdir'] + args.tem_input['b']
while (R0[0]<data.box[1, 0] and R0[1]<data.box[1, 1] and R0[2]<data.box[1, 2]):
t = t + dt
R0 = t*args.tem_input['vnorm']*args.tem_input['vdir'] + args.tem_input['b']
tmax = t-dt
t = np.arange(tmin, tmax+dt, dt)
r_elec = np.zeros((t.shape[0], 3), dtype=float)
r_elec[:, 0] = t*args.tem_input['vnorm']*args.tem_input['vdir'][0] + args.tem_input['b'][0]
r_elec[:, 1] = t*args.tem_input['vnorm']*args.tem_input['vdir'][1] + args.tem_input['b'][1]
r_elec[:, 2] = t*args.tem_input['vnorm']*args.tem_input['vdir'][2] + args.tem_input['b'][2]
return r_elec, t
#
#
#
def calc_force(rp, t, E_re, E_im):
"""
Calculate the force acting on the electron trajectory F = q.E,
where E is the induced field
"""
F = np.zeros(rp.shape, dtype=complex)
for i in range(t.shape[0]):
ind = pts_coord_to_Array_ind(E_re, rp[i, :])
F[i, :] = E_re.Array[ind[0], ind[1], ind[2], :] + complex(0.0, 1.0)*E_im.Array[ind[0], ind[1], ind[2], :]
return F
def calc_scatt_angle(args, V_re, V_im):
"""
Calculate the scaterring angle of the electron using the classical formulation
work only if the electron is going outside the sphere
"""
R_min = np.sqrt(np.dot(args.tem_input['b'], args.tem_input['b']))
theta = complex(0.0, 0.0)
for i in range(V_re.Array.shape[0]):
for j in range(V_re.Array.shape[1]):
for k in range(V_re.Array.shape[2]):
r = Array_ind_to_pts_coord(V_re, np.array([i, j, k]))
rnorm = np.sqrt(np.dot(r, r))
if rnorm > R_min:
y2 = (1-complex(V_re.Array[i, j, k], V_im.Array[i, j, k])/(0.5*args.tem_input['vnorm']**2))*rnorm**2
theta = theta + 1/(rnorm*np.sqrt(y2 - np.dot(args.tem_input['b'], args.tem_input['b'])))
return theta*V_re.dr[0]*V_re.dr[1]*V_re.dr[2]
#
#
#
def Force_electron(freq, vnorm, vdir, b, dt = 0.1, fname = 'tddft_tem_output.hdf5'):
"""
Calculate the force acting on the electron along its trajectory for a given frequency
Input parameters:
-----------------
freq (float): frequency at which calculate the force in eV
vnorm (float): electron velocity in au
vdir (float array (size=3)): velocity direction vector
b (float array (size=3)): impact parameter vector (au)
dt (float, default: 0.1): time step
fname (string, default: tddft_tem_output.hdf5): name of the input data, must be hdf5 file
"""
data = read_data()
data.args.quantity = "efield"
data.args.tem_iter = 'tem'
data.args.ReIm = 're'
data.args.plot_freq = freq
data.args.tem_input['vnorm'] = vnorm
data.args.tem_input['vdir'] = vdir
data.args.tem_input['b'] = b
E_re = data.Read(YFname=fname)
data.args.ReIm = 'im'
E_im = data.Read(YFname=fname)
r_elec, t = trajectory_elec(data.args, E_re, dt)
F_elec = calc_force(r_elec, t, E_re, E_im)
return F_elec, r_elec, t
#
#
#
def scattering_angle(args, fname = 'tddft_tem_output.hdf5'):
"""
Calculate the force acting on the electron along its trajectory for a given frequency
Input parameters:
-----------------
args
theta = pi - 2b\int_{R}^{+\infty}\frac{dr}{r}\frac{1}{\sqrt{y**2 - b**2}}
with
y(r) = r**2(1-\frac{V(r)}{E})
"""
data = read_data()
data.args = args
data.args.quantity = "potential"
data.args.tem_iter = 'tem'
data.args.ReIm = 're'
V_re = data.Read(YFname=fname)
data.args.ReIm = 'im'
V_im = data.Read(YFname=fname)
b = np.sqrt(np.dot(args.tem_input['b'], args.tem_input['b']))
theta = calc_scatt_angle(data.args, V_re, V_im)
return np.pi - 2*b*theta
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#initialisation of the pyiter package
from __future__ import division
import sys
from ase.cluster.icosahedron import Icosahedron
from ase.cluster.decahedron import Decahedron
from ase.cluster.octahedron import Octahedron
from ase.calculators.emt import EMT
from ase.optimize import BFGS
#
#
#
def constr_opt_cluster(symbol1, csym, nsh1, optimize=True, lcenter=True, lc=None):
if csym=="icosahedron" : cl = Icosahedron(symbol1, nsh1, latticeconstant=lc)
elif csym=="decahedron" : cl = Decahedron(symbol1, nsh1, nsh1, 0)
elif csym=="octahedron" : cl = Octahedron(symbol1, nsh1)
else :
print("unknown symmetry...")
sys.exit(11)
if lcenter: cl.positions = cl.positions - cl.get_center_of_mass()
cl.set_calculator(EMT())
E_start = cl.get_potential_energy()
if optimize :
print('Optimization starts at E_start', E_start)
dyn = BFGS(cl, trajectory="trajectory.traj", logfile="bfgs.log")
dyn.run(fmax=0.05)
E_finish = cl.get_potential_energy()
natoms = len(cl)
print('E_start/natoms, E_finish/natoms', E_start / natoms, E_finish / natoms)
return cl
#
#
#
def subst_atoms_cluster(cl_in, nsh2=1, symbol2="") :
from copy import deepcopy
cl = deepcopy(cl_in)
if symbol2=="" :
del cl [[atom.index for atom in cl if atom.tag<nsh2]]
else:
tags = cl.get_tags()
syms = cl.get_chemical_symbols()
for ia in range(len(cl)) :
if tags[ia]>nsh2 :
#print('tags[{0}]<{1}'.format(ia, nsh2))
syms[ia] = symbol2
cl.set_chemical_symbols(syms)
return cl
#
#
#
def constr_shell(symbol1="Ag", symbol2="", sym="icosahedron", nsh=2, nlayers=1, optimize=False, lc=None):
print('symbol2 = ', symbol2)
cl = constr_opt_cluster(symbol1, sym, nsh, optimize, True, lc)
cl1 = subst_atoms_cluster(cl, nsh - nlayers + 1, symbol2)
return cl1
#
#
#
def write_tags(fname, cl):
f=open(fname+'.tags', 'w')
f.write("{0}\n".format(len(cl)))
f.write(fname + ': atom -> tag\n')
for atom in cl:
f.write("{0} {1}\n".format(atom.index+1, atom.tag))
f.close()
def write_siesta_geo_ghost(atm, fname='geo'):
from ase.data import chemical_symbols
f=open(fname+'.fdf', 'w')
f.write('NumberOfAtoms {0}\n'.format(atm.get_number_of_atoms()))
f.write('NumberOfSpecies {0}\n'.format(atm.get_number_species()))
f.write('\n')
f.write('%block ChemicalSpeciesLabel\n')
sp_label = 1
sp = {}
nb = []
for i, atm_nb in enumerate(atm.get_atomic_numbers()):
if atm_nb not in nb:
f.write(' {0} {1} '.format(sp_label, atm_nb) + chemical_symbols[atm_nb]+'\n')
nb.append(atm_nb)
sp[chemical_symbols[atm_nb]] = sp_label
sp_label = sp_label + 1
if sp_label != atm.get_number_species()+1:
print('sp_label = {0} != nb_species = {1}'.format(sp_label, atm.get_number_species()))
raise ValueError('Something went wrong!!')
f.write('%endblock ChemicalSpeciesLabel\n')
f.write('\n')
f.write('AtomicCoordinatesFormat Ang\n')
f.write('%block AtomicCoordinatesAndAtomicSpecies\n')
ia = 0
for a in atm :
ia = ia + 1
f.write(" {0:.6f} {1:.6f} {2:.6f} {3} {4} ".format(a.x, a.y, a.z, sp[a.symbol], ia) + a.symbol+'\n')
f.write('%endblock AtomicCoordinatesAndAtomicSpecies\n')
print(fname+'.fdf')
f.close()
def write_siesta_geo(cl, fname='geo', ghost=None, ghost_atm_nb=-1):
from ase.data import chemical_symbols
f=open(fname+'.fdf', 'w')
f.write('NumberOfAtoms {0}\n'.format(len(cl)))
f.write('NumberOfSpecies {0}\n'.format(len(set(cl.numbers))))
f.write('%block ChemicalSpeciesLabel\n')
species_label={}
for i, nb in enumerate(set(cl.numbers)):
if chemical_symbols[nb] == 'X' and ghost is not None:
species_label[ghost] = i+1
f.write(' {0} {1} '.format(i+1, ghost_atm_nb)+ ghost+'\n')
elif chemical_symbols[nb] == 'X' and ghost is None:
raise ValueError('indicate the ghost specie')
else:
species_label[chemical_symbols[nb]] = i+1
f.write(' {0} {1} '.format(i+1, nb)+ chemical_symbols[nb]+'\n')
f.write('%endblock ChemicalSpeciesLabel\n')
f.write('AtomicCoordinatesFormat Ang\n')
f.write('%block AtomicCoordinatesAndAtomicSpecies\n')
ia = 0
for a in cl :
ia = ia + 1
if a.symbol == 'X' and ghost is not None:
f.write(" {0:.6f} {1:.6f} {2:.6f} {3} {4} ".format(a.x, a.y, a.z, species_label[ghost], ia) + ghost+'\n')
elif a.symbol == 'X' and ghost is None:
raise ValueError('indicate the ghost specie')
else:
f.write(" {0:.6f} {1:.6f} {2:.6f} {3} {4} ".format(a.x, a.y, a.z, species_label[a.symbol], ia) + a.symbol+'\n')
f.write('%endblock AtomicCoordinatesAndAtomicSpecies\n')
print(fname+'.fdf')
f.close()
return cl
#
#
#
def gen_siesta_geom_def(symbol1="Ag", symbol2="", sym="icosahedron", nsh=2, nlayers=1, optimize=False, lc=None, fname="geom"):
from ase.data import atomic_numbers
cl = constr_shell(symbol1, symbol2, sym, nsh, nlayers, optimize, lc)
f=open(fname+'.fdf', 'w')
f.write('NumberOfAtoms {0}\n'.format(len(cl)))
f.write('NumberOfSpecies {0}\n'.format(1))
f.write('%block ChemicalSpeciesLabel\n')
f.write(" {0} {1} ".format(1, atomic_numbers[symbol1])+ symbol1+"\n")
f.write('%endblock ChemicalSpeciesLabel\n')
f.write("\n")
f.write('AtomicCoordinatesFormat Ang\n')
f.write('%block AtomicCoordinatesAndAtomicSpecies\n')
ia = 0
for a in cl :
ia = ia + 1
sp = 1
f.write("{0:.6f} {1:.6f} {2:.6f} {3} {4} ".format(a.x, a.y, a.z, sp, ia) + a.symbol + "\n")
f.write('%endblock AtomicCoordinatesAndAtomicSpecies\n')
print(fname+'.fdf')
f.close()
return cl
#
#
#
def gen_siesta_geom_flt1(symbol1="Ag", symbol2="", sym="icosahedron", nsh=2, nlayers=1, optimize=False, lc=None, fname="geom"):
from ase.data import atomic_numbers
cl = constr_shell(symbol1, symbol2, sym, nsh+1, nlayers+1, optimize, lc)
f=open(fname+'.fdf', 'w')
f.write('NumberOfAtoms {0}\n'.format(len(cl)))
f.write('NumberOfSpecies {0}\n'.format(2))
f.write('%block ChemicalSpeciesLabel\n')
f.write(" {0} {1} ".format(1, atomic_numbers[symbol1]) + symbol1+"\n")
f.write(" {0} {1} ".format(2, -atomic_numbers[symbol1]) + symbol1+"s\n")
f.write('%endblock ChemicalSpeciesLabel\n')
f.write("\n")
f.write('AtomicCoordinatesFormat Ang\n')
f.write('%block AtomicCoordinatesAndAtomicSpecies\n')
ia = 0
for a in cl :
ia = ia + 1
sp = 1
if a.tag==nsh+1 : sp = 2
f.write("{0:.6f} {1:.6f} {2:.6f} {3} {4} ".format(a.x, a.y, a.z, sp, ia) + a.symbol + "\n")
f.write('%endblock AtomicCoordinatesAndAtomicSpecies\n')
print(fname+'.fdf')
f.close()
return cl
#
#
#
def gen_siesta_geom_flt2_acas(symbol1="Ag", symbol2="", sym="icosahedron", nsh=2, nlayers=1, optimize=False, lc=None, fname="geom"):
from ase.data import atomic_numbers
cl = constr_shell(symbol1, symbol2, sym, nsh+1, nlayers+2, optimize, lc)
f=open(fname+'.fdf', 'w')
f.write('NumberOfAtoms {0}\n'.format(len(cl)))
f.write('NumberOfSpecies {0}\n'.format(2))
f.write('%block ChemicalSpeciesLabel\n')
f.write("{0} {1} ".format(1 , atomic_numbers[symbol1]) + symbol1+"\n")
f.write("{0} {1} ".format(2 , -atomic_numbers[symbol1]) + symbol1+"s\n")
f.write('%endblock ChemicalSpeciesLabel\n')
f.write("\n")
f.write('AtomicCoordinatesFormat Ang\n')
f.write('%block AtomicCoordinatesAndAtomicSpecies\n')
ia = 0
for a in cl :
ia = ia + 1
sp = 1
if a.tag==nsh+1 : sp = 2
if a.tag==nsh-nlayers : sp = 2
f.write("{0:.6f} {1:.6f} {2:.6f} {3} {4} ".format(a.x, a.y, a.z, sp, ia) + a.symbol + "\n")
f.write('%endblock AtomicCoordinatesAndAtomicSpecies\n')
print(fname+'.fdf')
f.close()
return cl
#
#
#
def gen_siesta_geom_ico_ghost(symbol1, symbol2="", sym="icosahedron", nsh=2, nlayers=1, optimize=False, lc=None, fname="geom"):
from ase.data import atomic_numbers
if symbol2 == 'ghost':
cl = constr_shell("Ag", "X", sym, nsh+1, nlayers+2, optimize, lc)
else:
cl = constr_shell(symbol1, symbol2, sym, nsh+1, nlayers+2, optimize, lc)
f=open(fname+'.fdf', 'w')
f.write('NumberOfAtoms {0}\n'.format(len(cl)))
f.write('NumberOfSpecies {0}\n'.format(len(set(cl.numbers))))
f.write('%block ChemicalSpeciesLabel\n')
if symbol2 == 'ghost':
species_label={symbol1: [atomic_numbers[symbol1], 1],
symbol1+"s": [-atomic_numbers[symbol1], 2]}
for k, val in species_label.items():
f.write(' {0} {1} '.format(val[1], val[0])+ k+'\n')
else:
species_label={symbol1: [atomic_numbers[symbol1], 1],
symbol2: [atomic_numbers[symbol2], 2]}
for k, val in species_label.items():
f.write(' {0} {1} '.format(val[1], val[0])+ k+'\n')
f.write('%endblock ChemicalSpeciesLabel\n')
f.write('ZM.UnitsLength Ang\n')
f.write('%block Zmatrix\n')
f.write('cartesian\n')
ia = 0
for a in cl :
ia = ia + 1
ox = 1; oy = 1; oz = 1;
if a.tag==nsh+1 or a.tag==nsh-nlayers :
ox = 1; oy = 1; oz = 1
#print(a.symbol)
if a.symbol == 'X':
symbol = symbol1+"s"
else:
symbol = a.symbol
f.write("{0} {1:.5f} {2:.5f} {3:.5f} {4} {5} {6} {7} ".format(species_label[symbol][1],
a.x, a.y, a.z, ox, oy, oz, ia)+symbol+'\n')
f.write('%endblock Zmatrix')
print(fname+'.fdf')
f.close()
return cl
#
#
#
def gen_siesta_geom_flt2_zmat(symbol1, symbol2="", sym="icosahedron", nsh=2, nlayers=1, optimize=False, lc=None, fname="geom"):
from ase.data import chemical_symbols
cl = constr_shell(symbol1, symbol2, sym, nsh+1, nlayers+2, optimize, lc)
f=open(fname+'.fdf', 'w')
f.write('NumberOfAtoms {0}\n'.format(len(cl)))
f.write('NumberOfSpecies {0}\n'.format(len(set(cl.numbers))))
f.write('%block ChemicalSpeciesLabel\n')
species_label={}
for i, nb in enumerate(set(cl.numbers)):
species_label[chemical_symbols[nb]] = i+1
f.write(' {0} {1} '.format(i+1, nb)+ chemical_symbols[nb]+'\n')
f.write('%endblock ChemicalSpeciesLabel\n')
f.write('ZM.UnitsLength Ang\n')
f.write('%block Zmatrix\n')
f.write('cartesian\n')
ia = 0
for a in cl :
ia = ia + 1
ox = 1; oy = 1; oz = 1;
if a.tag==nsh+1 or a.tag==nsh-nlayers :
ox = 1; oy = 1; oz = 1
f.write("{0} {1:.5f} {2:.5f} {3:.5f} {4} {5} {6} {7} ".format(species_label[a.symbol], a.x, a.y, a.z, ox, oy, oz, ia)+\
a.symbol+'\n')
f.write('%endblock Zmatrix')
print(fname+'.fdf')
f.close()
return cl
#
# Main prog
#
#symbol1 = "Ag"
#symbol2 = ""
#optimize = False
#lcenter = True
#sym = "icosahedron"
#lc = 5.07
#
#for nsh in [2, 3, 4, 5, 6, 7] :
# ll = list(set([nsh]))
# for nlayers in ll :
# if(nlayers>nsh): continue
# fname = sym[0]+"_sz"+str(nsh)+"_la"+str(nlayers)+"_"+symbol1+symbol2+\
# "_flt2_zmat"
# cl = gen_siesta_geom_flt2_zmat(symbol1, symbol2, sym, nsh, nlayers, \
# optimize, lc=lc, fname='geometry.siesta')
#
#
#
# ase.io.write(fname+'.xyz', cl)
# print fname+'.xyz'
# write_tags(fname, cl)
# print fname+'.tags'
# folder = 'Na{0}_lc{1:.2f}'.format(cl.get_positions().shape[0], lc)
# run('mkdir '+folder)
# run('mv '+fname+'.* '+folder)