#!/usr/bin/env python3.9

'''

#=======================================================================

#Simple script to calculate the Madelung potential (MP)

#

#Author : Asif Iqbal -> @AIB_EM

#Created on : 14/05/2020

#

#USAGE : python3 argv.sys[0] It will read the CONTCAR automatically.

#NB: NO warranty guarranteed whatsoever or even implied. The script reads

#CONTCAR file and compute the Madelung potential. It is better to read the final

#position after relaxation and associated charges.

#---

#This is a simple script which assumes the charges on the atom with atomic charges.

#In real case it should be approximated with Mulliken charge or Bader charge analysis.

#The charges can be appended to the 4th column in the CONTCAR file and the script

#can be modified to read the charges on each atom and calculate the MP.

#=======================================================================

'''

import sys, os

import numpy as np

import warnings

# Constants

E = 1.602176634e-19 # Elementary charge (Coulombs)

P = 1.11265e-10 # 4πϵo is 1.11265x10-10 C**2/(Jm).

A2M = 1e-10 # Angstrom to meter

J2eV = 6.242e+18

contcar_file = 'POSCAR'

if not os.path.exists(contcar_file):

print('ERROR: CONTCAR does not exist')

sys.exit(1)

print('Reading CONTCAR ...\n')

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

lines_contcar = f.readlines()

for line in lines_contcar:

if line.strip().lower() in ["direct", "d"]:

warnings.warn('The file should be in Cartesian coordinates. Convert to Cartesian before proceeding.')

sys.exit(1)

#-------------------------------------------------------------------

# Defining atomic numbers

#-------------------------------------------------------------------

# THE CHEMICAL SYMBOL OF ELEMENTS IN THE PERIODIC TABLE,

# EXTRACTED FROM VESTA CONFIGURATION FILE.

atomic_name = [

"Np", "Pu", "Am", "XX"]

atomic_number = [

92, 93, 94, 95, 96]

atomic_charge = [

+3, 0, 0, 0, 0]

atomic_mass = {

'Cr': 51.996, 'Ca': 40.08, 'Cu': 63.546, 'In': 114.82}

#-------------------------------------------------------------------

dict, at_nu = {}, {}

rx, ry, rz = [], [], []

qq, r, charge = [], [], []

#============== READING ATOM TYPES AND # OF ELEMENTS

atoms_types, atom_list = lines_contcar[5].split(), lines_contcar[6].split()

output = f"{'ATOMS TYPE':30s}: {list(zip(atoms_types, atom_list))}"

print(output)

#==============

# SUMMING ATOM NUMBERS

#==============

sum_atoms = atom_list

for i in range( len(atoms_types) ):

dict[i] = atoms_types[i], sum_atoms[i]

sum_atoms = sum(list(map(int, sum_atoms)))

print(f"{'TOTAL # OF ATOMS':30s}: {sum_atoms}")

#==============

# FINDING ATOMIC NUMBER FROM A LIST &

# MULTIPLYING THE LIST WITH ATOMIC NUMBERS

# ACCORDING TO THE POSCAR LIST FORMAT

#==============

for j in range(len(atoms_types)):

for k in atomic_name:

if dict[j][0] == k:

element = atomic_name.index(k)

at_nu[j] = dict[j][0], atomic_number[element]

#print(dict[j][0],"-->", atomic_number[element])

#print (f'Charge Dictionary -> {at_nu}')

for i in range( len(atoms_types) ):

charge.extend(at_nu[i][1] for _ in range( int(atom_list[i]) ))

print(f"{'Charge Multiplicity':30s}: {charge}")

#============== INDEXING POSITIONS

for i in lines_contcar:

if any(s in i for s in ["Cartesian", "C", "cartesian"]):

lp = lines_contcar.index(i)

coord_q = np.array([list(map(float, lines_contcar[lp+1+i].split())) for i in range(sum_atoms)])

rx, ry, rz, qq = coord_q.T

print(f"{'='*50:<30}")

print(f"CHARGES IN FILE {qq}\n")

temp = 0.0

for i in range(sum_atoms):

for j in range(sum_atoms):

if (j > i): # AVOIDING DOUBLE COUNTING

r = np.sqrt( (rx[i]-rx[j])**2 + (ry[i]-ry[j])**2 + (rz[i]-rz[j])**2 ) * A2M

temp += ( qq[j]/r )

break

Medulung = temp * (E/P)

print(f"Medulung Potential is: {Medulung:6.3f} J" )