2: (2.02290, 0.00000),

3: (0.00000, 1.12780),

4: (2.02290, 1.12780),

5: (0.00000, 2.25560),

6: (2.02290, 2.25560),

7: (0.00000, 3.38340),

8: (2.02290, 3.38340),

9: (0.00000, 4.51120),

2: [],

3: [],

4: [],

5: [],

6: [],

7: [],

8: [],

9: [],

2: [],

3: [],

4: [],

5: [],

6: [],

7: [],

8: [],

9: [],

2: [],

3: [],

4: [],

5: [],

6: [],

7: [],

8: [],

9: [],

1: '#e34545',

2: '#de27e5',

3: '#5d4fe5',

4: '#4fb9e5',

5: '#10be3d',

6: '#3f5d36',

7: '#7d4804',

8: '#ff783c',

9: '#c70707',

10: '#560606'

"""

poscar = "POSCAR.1"

with open(poscar, 'r') as f:

# Skipping the comment

f.readline()

acell = float(f.readline())

cell = []

for i in range(0, 3):

cell.append(f.readline().split())

cell = np.array(cell) # Cell will be an array of strings

cell = cell.astype(np.float) # Transform cell into an array of floats

# read kind of atoms

at_kinds = f.readline().split()

# read number of each kind of atoms

no_at_each_kind = f.readline().split()

total_atoms = 0

for i in no_at_each_kind:

total_atoms = total_atoms + int(i)

f.readline() # Skip 'selective dynamics' line

coord_type = f.readline()

# If coordinates are in direct coordinates, transform to cartesians

transform = False

if coord_type[0].lower() == "d":

transform = True

idx = 1

coordinates = []

constraints = []

while idx <= total_atoms:

line = f.readline()

if line != "\n":

coordinates_i = line.split()[0:3]; constraints_i = line.split()[3:7]

if transform:

coordinates_i = transform_coordiantes(cell, coordinates_i, acell)

coordinates.append([idx, coordinates_i])

constraints.append([idx, constraints_i])

idx += 1

atom_index = 1

z = coordinates[0][1][2]

for atom in coordinates:

if kind == 'tip':

# if atom[1][2] > float(z):

if atom[1][2] > z:

atom_index = atom[0]

z = atom[1][2]

else:

# if atom[1][2] < float(z):

if atom[1][2] < z:

atom_index = atom[0]

z = atom[1][2]

# print (kind + " index: ", atom_index)

# ------ | |--> z_tip |

# \ / |l_tip | |

# \ / | | |

# \/ | | |

# | | h | z = (z_tip - z_surf) - l

# | | | '------.-------'

# z | | | h

# | | |

# __|______ | |

# ///////// |l_surf |

# ///////// | |_,> z_surf

min_z = 1000.0

max_z = -1000.0

for d in coordinates:

z = float(d[1][2])

if z > max_z:

max_z = z

if z < min_z:

min_z = z

distance = abs(max_z - min_z)

# for i in range (0, len(coordinates)):

# idx = coordinates[i][0]

# if idx == at_1:

# coord_at_1 = coordinates[i][1]

# if idx == at_2:

# coord_at_2 = coordinates[i][1]

# # min_z = coordinates[0][1]

# # max_z = coordinates[0][1]

# # for i in range (0, len(coordinates)):

# distance = abs(float(coord_at_1[2]) - float(coord_at_2[2]))

distance = distance - tip_length - surface_length

""" Transform coordinates from Direct to Cartesian """

# Convert the values to numpy arrays

cell = np.array(cell)

coordinates = np.array(coordinates)

coordinates = coordinates.astype(np.float)

new_coordinates = []

if coordinates.ndim == 1:

x = coordinates[0] * cell[0, :] * acell

y = coordinates[1] * cell[1, :] * acell

z = coordinates[2] * cell[2, :] * acell

new_coordinates = x + y + z

return new_coordinates

else:

# TODO: This block does not work.

# If coordinates is an array of coordinates of several atoms, the

# function does not work

for i in range(0, len(coordinates)):

x = coordinates[i, 0] * cell[0, :] * acell

y = coordinates[i, 1] * cell[1, :] * acell

z = coordinates[i, 2] * cell[2, :] * acell

coord_i = x + y + z

new_coordinates.append(coord_i)

""" Reads the POSCAR (for VASP 5.x)"""

with open(poscar, 'r') as f:

# Skipping the comment

f.readline()

acell = float(f.readline())

cell = []

for i in range(0, 3):

cell.append(f.readline().split())

cell = np.array(cell) # Cell will be an array of strings

cell = cell.astype(np.float) # Transform cell into an array of floats

# read kind of atoms

at_kinds = f.readline().split()

# read number of each kind of atoms

no_at_each_kind = f.readline().split()

total_atoms = 0

for i in no_at_each_kind:

total_atoms = total_atoms + int(i)

f.readline()

coord_type = f.readline()

# If coordinates are in direct coordinates, transform to cartesians

transform = False

if coord_type[0].lower() == "d":

transform = True

# print "coordinates in file", poscar, "are Direct. Will be transformed to cartesian"

idx = 1

coordinates = []

constraints = []

while idx <= total_atoms:

line = f.readline()

if line != "\n":

coordinates_i = line.split()[0: 3]; constraints_i = line.split()[3: 7]

if transform:

coordinates_i = transform_coordiantes(cell, coordinates_i, acell)

coordinates.append([idx, coordinates_i])

constraints.append([idx, constraints_i])

idx += 1

global f_interp

if method == "poly":

# Fit Distance-Energy curve to a polynomy of rank <pol_rank>

# Overwriting pol_rank value (shall I keep it in this way?)

pol_rank = 12

p = np.poly1d(np.polyfit(distances, energies, pol_rank))

# Compute the -derivative of the polynomy (the force)

pp = -np.polyder(p)

force = pp

f_interp = force

# return p, force

if method == "pol_interp4":

force = []

for i in range(4, len(distances)-4):

section_x = distances[i-4:i+5]

section_y = energies[i-4:i+5]

# Overwriting pol_rank value (shall I keep it in this way?)

pol_rank = 5

# Fit section of Distance-Energy curve to a polynomy of rank <pol_rank>

p = np.poly1d(np.polyfit(section_x, section_y, pol_rank))

# Compute the -derivative of the polynomy (the force)

pp = -np.polyder(p)

f_i = np.polyval(pp, distances[i])

force.append(f_i)

f_interp = force

# return force

if method == "pol_interp3":

force = []

for i in range(3, len(distances)-3):

section_x = distances[i-3:i+4]

section_y = energies[i-3:i+4]

# Overwriting pol_rank value (shall I keep it in this way?)

pol_rank = 5

# Fit section of Distance-Energy curve to a polynomy of rank <pol_rank>

p = np.poly1d(np.polyfit(section_x, section_y, pol_rank))

# Compute the -derivative of the polynomy (the force)

pp = -np.polyder(p)

f_i = np.polyval(pp, distances[i])

force.append(f_i)

f_interp = force

# return force

if method == "pol_interp2":

force = []

for i in range(2, len(distances)-2):

section_x = distances[i-2:i+3]

section_y = energies[i-2:i+3]

# Overwriting pol_rank value (shall I keep it in this way?)

pol_rank = 5

# Fit section of Distance-Energy curve to a polynomy of rank <pol_rank>

p = np.poly1d(np.polyfit(section_x, section_y, pol_rank))

# Compute the -derivative of the polynomy (the force)

pp = -np.polyder(p)

f_i = np.polyval(pp, distances[i])

force.append(f_i)

f_interp = force

# return force

if method == "pol_interp":

force = []

for i in range(1, len(distances)-1):

section_x = distances[i-1:i+2]

section_y = energies[i-1:i+2]

# Overwriting pol_rank value (shall I keep it in this way?)

pol_rank = 5

# Fit section of Distance-Energy curve to a polynomy of rank <pol_rank>

p = np.poly1d(np.polyfit(section_x, section_y, pol_rank))

# Compute the -derivative of the polynomy (the force)

pp = -np.polyder(p)

f_i = np.polyval(pp, distances[i])

force.append(f_i)

f_interp = force

# return force

if method == "num1":

# Compute the forces via numerical derivative

dE = []

for i in range(0, len(distances)-1):

dE_i = -(energies[i]-energies[i+1])/(distances[i]-distances[i+1])

dE.append(dE_i)

dE_fit = np.poly1d(np.polyfit(distances[:-1], dE, pol_rank))

force = dE

f_interp = force

# return force

if method == "num2":

# Using (f(x-h) - f(x+h)) / 2h

dE_2 = []

for i in range(1, len(distances)-1):

dE_i = -(energies[i-1]-energies[i+1])/(distances[i-1]-distances[i+1])

dE_2.append(dE_i)

dE_2_fit = np.poly1d(np.polyfit(distances[1:-1], dE_2, pol_rank))

force = dE_2

f_interp = force

# return force

if method == "spline":

global force_spl

global e_fit_spl

e_fit_spl = None

force_spl = None

# Fit the curve E(z) to a spline

e_fit_spl = UnivariateSpline(distances, energies, s=spl_smoothing)

# Compute the first derivative of the spline

force_spl = e_fit_spl.derivative(1) # Remember to use -force_spl later

"""Read OUTCAR type file.

if isinstance(filename, str):

f = open(filename)

else: # Assume it's a file-like object

f = filename

data = f.readlines()

natoms = 0

images = []

# atoms = Atoms(pbc = True)

energy = 0

atoms = []

forces = []

species = []

symbols = []

species_num = []

for n, line in enumerate(data):

if 'VRHFIN' in line:

temp = line.split('=')

species.append(temp[1][0:2].strip(':'))

if 'ions per type' in line:

temp = line.split()

for ispecies in range(len(species)):

species_num += [int(temp[ispecies+4])]

natoms += species_num[-1]

for iatom in range(species_num[-1]): symbols += [species[ispecies]]

if 'direct lattice vectors' in line:

cell = []

for i in range(3):

temp = data[n+1+i].split()

cell += [[float(temp[0]), float(temp[1]), float(temp[2])]]

if 'FREE ENERGIE OF THE ION-ELECTRON SYSTEM' in line:

energy = float(data[n+2].split()[4])

if 'POSITION ' in line:

# Restart here the forces and atoms list in order to take just the last geometry

forces = []

atoms = []

for iatom in range(natoms):

temp = data[n+2+iatom].split()

atoms += [ (symbols[iatom], [float(temp[0]), float(temp[1]), float(temp[2])])]

forces += [[float(temp[3]), float(temp[4]), float(temp[5])]]

images = [(atoms, forces)]

# Uncomment to store all the iterations

# images += [(atoms, forces)]

# TODO: Check that the OUTCAR has converged

write_geom_and_forces(natoms, atoms, forces)

jmol_script = 'jmolscript: vectors on; vectors 2; set vectorscale 2; set percentVdwAtom 30; set bondradiusmilliangstroms 120;'

fout = open('temp.xyz', 'a')

fout.write(str(natoms) + "\n")

fout.write(jmol_script + "\n")

for i in range (0, len(atoms)):

# Sorry for the long line

fout.write('%s %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f \n' % ( atoms[i][0], atoms[i][1][0], atoms[i][1][1], atoms[i][1][2], forces[i][0], forces[i][1], forces[i][2]))

atoms, forces = read_outcar (outcar)

vertical_force = 0

for i in range(0, len(forces)):

if constraints[i][1][2].lower() == "f":

# print atoms[1] , forces[i], constraints[i][1]

# print 'forces: ', forces[i][2], constraints[i][1][2]

vertical_force += forces[i][2]

'''

'''Set a E fitting interval between z_min and z_max'''

# Select points between height (z) interval

global dist_interval

global ener_interval

global outcar_interval

# Delete previously set interval

del dist_interval[:]

del ener_interval[:]

del outcar_interval[:]

dist = []

ener = []

outc = []

for i, z in enumerate(distances):

if z >= z_min and z <= z_max:

# dist.append(distances[i])

# ener.append(energies[i])

# outc.append(outcar_files[i])

# distance

dist.append(spectroscopy[i][0])

# energy

ener.append(spectroscopy[i][1])

# Outcar file

outc.append(spectroscopy[i][2])

dist_interval = dist

ener_interval = ener

'''Plot outcar labels'''

# Add labels to energy points

idx_label = 0

for i in range(len(ener_interval)):

if (idx_label % label_interval == 0 ):

ax1.text(dist_interval[i]+0.2, ener_interval[i], outcar_interval[i].replace("OUTCAR.",""),fontsize=8)

ax2.text(dist_interval[i]+0.2, -force_spl(dist_interval[i]), outcar_interval[i].replace("OUTCAR.",""), fontsize=8)

ax3.text(dist_interval[i]+0.2, ener_interval[i], outcar_interval[i].replace("OUTCAR.",""), fontsize=8)

''' Replot graphs'''

# Clear previous data but not axes

ax1.cla()

ax2.cla()

ax3.cla()

# Padding

plt.tight_layout(pad=1.1)

ax1.margins(0.05, 0.2)

ax2.margins(0.05, 0.2)

ax3.margins(0.05, 0.2)

#plot_name = prefix+"_E_"+str(grid_point)+".pdf"

ax1.plot(dist_interval, ener_interval, 'o', label="Energies", lw=1)

ax1.plot(dist_interval, e_fit_spl(dist_interval), '--', label="E fit", lw=1)

#plot_name = prefix+"_F_"+str(grid_point)+".pdf"

#ax2.plot(dist_interval, -force(dist_interval), 'k-', label="Force", lw=1)

# Computed via with spline

ax2.plot(dist_interval, -force_spl(dist_interval), 'k-', label="Force", lw=1)

#plot_name = prefix+"delta_w_"+str(grid_point)+".pdf"

ax3.plot (h,delta_w, 'm-', label = 'omega', lw=1)

#ax3.plot (h,delta_w, 'm-', lw=1.5)

# Energy

ax1.set_title('Energy', fontsize=10)

ax1.set_xlabel (r'z ($\AA$)')

ax1.set_ylabel ('Energy (ev)')

# Force

ax2.set_title('Force vs vertical distance', fontsize=10)

ax2.set_xlabel (r'z ($\AA$)')

ax2.set_ylabel (r'Force(ev/$\AA$)')

# Omega

ax3.set_title('Frequency shift', fontsize=10)

ax3.set_xlabel (r'h ($\AA$)')

ax3.set_ylabel (r'$\Delta \omega$ (Hz)')

#plot_name = prefix+"_allgraphs_"+str(grid_point)+".pdf"

if add_labels: plot_outcar_labels()

'''Calculate the integral of the forces'''

global w

global h

global delta_w

# Empty the values before computing the integrals

w = [] # w => (w_exp/w_0)**2

h = [] # h : z + and - the apmplitud

delta_w = [] # (w_exp - w_0)

w_vs_h = []

min_z = dist_interval[0]

max_z = dist_interval[-1]

# Integrate at every z value

for hh in dist_interval:

if hh-A < min_z: continue # Check that we are not integrating out of the fiting boundaries

if hh+A > max_z: continue

# NOTE: With this xx value, got good results

xx = np.linspace(hh-A, hh+A, int_steps)

# NOTE: With this xx value all integrals are 0

# xx = hh + A*np.sin(phi)

#f_sum = -force_spl(hh + np.sin(phi))* np.sin(phi)*(2*np.pi/int_steps)

# f_int_sum=sum(f_sum)

f = -force_spl(hh + A*np.sin(phi)) * np.sin(phi)

# # KKKK Temp file to write out the forces and read in into fortran code

# fla = open('fla.dat_'+str(hh), 'w')

# fla.write('# xx hh+sin(phi) f(h+Asin(phi))sin(phi) phi \n')

# for i, fi in enumerate(f):

# fla.write(str(xx[i]) + "\t" + str(hh + A*np.sin(phi[i])) + "\t" + str(fi) + "\t" + str(phi[i]) + '\n')

f_int = simps(f, phi)

# fla.write("# Integral: " + str(f_int))

# fla.close()

function = 1.0-1.0/(np.pi*k*A)*f_int

w_vs_h += [(hh, function)]

for i in range(len(w_vs_h)):

h.append(w_vs_h[i][0]); w.append(w_vs_h[i][1])

delta_w = []

for i, w_i in enumerate(w):

temp = np.sqrt(w_i)*w_0

''' Read autosaved file. Note ifile is an already read file as f.readlines()'''

for n, line in enumerate(ifile):

if 'grid_point' in line and not 'end' in line:

if int(grid_point) == int(line.split()[1]):

# range 4: Max number of options expected in this block

for i in range(4):

option = ifile[n+i]

print(option)

if 'end_grid_point' in option:

break

elif 'smoothing' in option:

global spl_smoothing

spl_smoothing = float(option.split()[1])

elif 'z_min' in option:

global z_min

z_min = float(option.split()[1])

elif 'z_max' in option:

global z_max

global spl_smoothing

global z_min

global z_max

print(grid_point)

ofile.write("grid_point " + str(grid_point) + '\n')

ofile.write(" z_min " + str(z_min) + '\n')

ofile.write(" z_max " + str(z_max) + '\n')

ofile.write(" smoothing " + str(spl_smoothing) + '\n')

''' Read input file. Default inp.afm'''

if isinstance(filename, str):

f = open(filename)

else: # Assume it's a file-like object

f = filename

data = f.readlines()

for n, line in enumerate(data):

if line.startswith('#'):

continue

if 'smoothing' in line:

global spl_smoothing

spl_smoothing = float(line.split()[1])

elif 'z_min' in line:

global z_min

z_min = float(line.split()[1])

elif 'z_max' in line:

global z_max

z_max = float(line.split()[1])

elif 'omega_exp' in line:

global omega_exp

omega_exp = float(line.split()[1])

elif 'tip_length' in line:

global tip_length

tip_length = float(line.split()[1])

elif 'surface_length' in line:

global surface_length

surface_length = float(line.split()[1])

elif 'datfile' in line:

global outdatfile

outdatfile = str(line.split()[1])

elif 'tip_name' in line:

global tip_name

# The tip_name will be up to the end of the line

tip_name = line[9::]

elif 'add_labels' in line:

global add_labels

add_labels = string_to_bool(line.split()[1])

elif 'label_interval' in line:

global label_interval

label_interval = int(line.split()[1])

elif 'prefix' in line:

global prefix

prefix = str(line.split()[1])

elif 'grid_points' in line:

global plot_grid_points

plot_grid_points = []

for i in line.split()[1::]:

plot_grid_points.append(int(i))

elif 'show_graphs' in line:

global show_graphs

show_graphs = string_to_bool(line.split()[1])

elif 'save_forces' in line:

global save_forces

save_forces = string_to_bool(line.split()[1])

elif 'save_graphs' in line:

global save_graphs

save_graphs = string_to_bool(line.split()[1])

elif 'include_retractions' in line:

global include_retractions

include_retractions = string_to_bool(line.split()[1])

elif 'only_forces' in line:

global only_forces

only_forces = string_to_bool(line.split()[1])

elif 'do_full_image' in line:

global do_full_image

do_full_image = string_to_bool(line.split()[1])

elif 'save_afmimage' in line:

global save_afmimage

save_afmimage = string_to_bool(line.split()[1])

elif 'show_afmimage' in line:

global show_afmimage

show_afmimage = string_to_bool(line.split()[1])

elif 'autosave' in line:

global autosave

autosave = string_to_bool(line.split()[1])

elif 'force_initial' in line:

global force_initial

force_initial = string_to_bool(line.split()[1])

elif 'overwrite' in line:

global overwrite

overwrite = string_to_bool(line.split()[1])

else:

global afmimage

global grid

# GP(x,y) coord, h vs delta_w

# | | |

global afmforces

global grid

# GP(x,y) coord, h vs forces

# | | |

f = np.array(curve)-omega

ffit = UnivariateSpline(h, f, s=0)

roots = ffit.roots()

if roots.size > 1:

#root = max(roots) # Take the Delta_w from the "atractive" part => Gives no contrast

root = min(roots) # Take the Delta_w from the "repulsive" part => Gives contrast

elif roots.size == 1:

root = roots

else:

print("Error! No Z found at this delta Omega. Please increase your Delta Omega and try again")

return