Difference in signal for the same molecule using particleAtoms (with PDB-file) and ParticleMap (with meshgrid) at q=0

[idawod]

Hi,

Comparing the resulting signal from the same molecule using two different methods, particleAtoms and particleMap, I do not get the same value for q=0 (If I don't scale with a factor). The input for particleAtoms is a PDB-file. Using the same PDB-file, I have with DFT calculated the electron density on a grid, and made it an input to Condor. All files to run on the cluster is attached.

PDB_vs_DENCHAR.zip

[mhantke]

Hi,

there is a test in tests.py that should check for that kind of consistency. Have you compared your code with the one in the test?

Cheers, Max

[FilipeMaia]

I think the issue is also about understanding what condor.ParticleMap does when used with map3d and electron_density. It creates a material with an average electron density as the argument given and from it calculates a complex index of refraction. It then uses the map3d to scale that complex index of refraction on each voxel.

[idawod]

I had a look at those tests, and tried them myself with my input. At the time I tested, it did not pass the test. But the input I used is a bit different than the map used in the tests.

my input looks like this:

particle = condor.ParticleMap(map3d=charge_density, electron_density=3.332e28*10, dx = DX, geometry="custom", material_type="custom", rotation_values=rotation, rotation_formalism="quaternion")

The charge_density input to map3d is my electron density from the simulations, so it is multiplied by 1/electron_density.

I integrated my meshgrid input (charge_density) using the gridspacing (dx) and np.trapz. The result gives a value very close to 92, which is the number of electrons in the system for tri-alanine.

From what I understand, both methods should be exactly the same at q=0. If you run the scripts, the averaged intensity as a function of q has a similar shape up to a certain q (but the scaling is off).

I'll attach a test I did with a larger molecule: Trp-cage (see figure). And also a figure where I have scaled the intensity such that the to the methods are the same at q=0:

[mhantke]

Are you sure that you correctly rescaled the output of your parse_cube script to electron densities (unit has to be 1/m3)? I get a discrepancy by a factor of almost exactly 5.

# Read .cube file
atoms, electron_density_bohr, gridspacing = parse_cube("neutral_250_gridpoints.cube")

# Count electrons from atoms
atom_to_z = {
    "D": 1,
    "C": 6,
    "N": 7,
    "O": 8,
} 
n_electrons_from_atoms = np.asarray([atom_to_z[atom.name] for atom in atoms]).sum()
print("Number of electrons (from atoms): n(atoms) = %i" % n_electrons_from_atoms)

# Count electrons from density
bohr = 0.529177249E-10
dx = gridspacing[0, 0] * bohr
dy = gridspacing[1, 1] * bohr
dz = gridspacing[2, 2] * bohr
dv = dx * dy * dz
electron_density = electron_density_bohr / bohr**3
n_electrons_from_density = electron_density.sum() * dv
print("Number of electrons (from density): n(density) = %f" % n_electrons_from_density)

print("n(density) / n(atoms) = %f" % (n_electrons_from_density/n_electrons_from_atoms))

Output:

Number of electrons (from atoms): n(atoms) = 124
Number of electrons (from density): n(density) = 620.843958
n(density) / n(atoms) = 5.006806

[idawod]

Hi,

the grid-spacing is outputed in bohr, and the mesh-grid in angstrom. If you multiply the "Number of electrons (from density): n(density) = 620.843958" by bohr**3, you will get 92.

However, I think I have realized what the issue is, and it is not a problem of Condor. The fact that you calculated the number of electrons to be 124, made me have a look at my DFT simulations again. The software used for the simulations uses pseudo-potentials for the atoms: so it does not explicitly calculate the core electron density (to speed up calculations). The density corresponding to the mesh-grid is for the valence electrons only. The total number of electrons in the molecule is 124, and without core electrons the number is 92, which results in a factor: 124/92 = 1.3478, which squared becomes approximately 1.81, corresponding to the factor needed to be multiplied at q=0 to get the same value for the two methods.

Since the output of the DFT simulations gave the total number of electrons: 92, and I got the same by integrating the mesh-grid, I did not think I was missing electrons. But, I have realized now that it does not take the core electron density into account when writing out the output data.

So I think we can probably close this issue. Thanks a lot for taking the time to help!

Cheers

[FilipeMaia]

@idawod Glad you figured it out!