Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

comer mann electron dens #38

Open
dermen opened this issue Oct 9, 2018 · 0 comments
Open

comer mann electron dens #38

dermen opened this issue Oct 9, 2018 · 0 comments

Comments

@dermen
Copy link
Collaborator

dermen commented Oct 9, 2018

It seems that the scatter.atomic_electrondens had a small typo:

phi[inds] = cromermann[i] * \

should be

phi[inds] += cromermann[i] * \

Other than that there was a discrepancy from the formula I found via quick wolfram, and the one used.. how does this look:

click here for wolfram alpha computation

I tried loading the referenced paper in the function documentation, but could not find it (minimal effort used)

I made a little test script :

import numpy as np
import pylab as plt
from thor import scatter
def single_atom_dens(atomic_num, grid_size=5., grid_spacing=0.05):
    """
    atomic_num, atomic number of element 
    grid_size, size of cubic density grid in Angstrom
    grid_spacing, spacing of cubic grid in Angstrom
    """
    Ngrid_pt = int( float(grid_size)/grid_spacing)

    vals = np.arange( -grid_size/2, grid_size/2, grid_spacing)
    
    x,y,z = np.meshgrid( vals,vals,vals, sparse=True)
    R = np.sqrt( x**2 + y**2 + z**2 )
    form_fact = scatter.atomic_electrondens( atomic_num, R.ravel() )

    return form_fact.reshape( (Ngrid_pt, Ngrid_pt, Ngrid_pt))


# get the density
some_dens = {n:single_atom_dens(n) for n in range( 1,26)}


# plot
fig,axs = plt.subplots( 5,5, figsize=(10,10))
atom = 1
for i in range (5):
    for j in range(5):
        axs[i][j].clear()
        axs[i][j].imshow( sum(some_dens[atom],0), vmax=50)
        axs[i][j].set_title( "Z=%d"%atom, fontsize=15)
        axs[i][j].axis('off')
        atom += 1
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant