Dplot output charge density, total number of electron incorrect
Just Got Here
8:42:16 AM PDT - Wed, Oct 16th 2013
Dear NwChem Users,

I want to calculate the total charge density in a cube file by mp2 method. While, after I have this done, I integrate all the charge density, but it is not the total number of electron in the system. Is there something wrong with my input file?

Here is my input file and sum charge script.

Yi Yao

mp2.nw

=

start nacl

geometry
 na 0 0 0

 cl 2.0 0 0
end

basis
 * library 6-311G**
end

mp2
 freeze core
end

task mp2

dplot
 TITLE CHARGE_DENSITY

 vectors nacl.movecs

  LimitXYZ

-7.5 7.5 100

-7.5 7.5 100

-7.5 7.5 100

 spin total

 gaussian

 output nacl.cube
end

task dplot


sum_charge.py

=

import sys,math

cube_file_name = sys.argv[1]

cube_file = file(cube_file_name)

init =[0,0,0]
npix =[0,0,0]
dpix =[0,0,0]
sumcharge = 0

          1. READ THE HEAD INFORMATION
cube_file.readline()
cube_file.readline()
line = cube_file.readline()
words = line.split()
natoms = int(words[0])
for i in range(3):
   init[i] = float(words[i+1])
for i in range(3):
   line = cube_file.readline()

   words = line.split()

   npix[i] = int(words[0])

   dpix[i] = float(words[i+1])
  1. print init, npix, dpix,
for i in range(natoms):
   cube_file.readline(),

for i in range(npix[0]):
   #print i

   for j in range(npix[1]):

       for k in range(int(math.ceil(npix[2]/6.0))):

           line = cube_file.readline()

           words = line.split()

           for word in words:

               sumcharge += float(word)
total_charge = sumcharge*dpix[0]*dpix[1]*dpix[2]
print total_charge
cube_file.close()

Forum Vet
2:03:39 PM PDT - Wed, Oct 16th 2013
Quote:Yiy Oct 16th 7:42 am
Dear NwChem Users,

I want to calculate the total charge density in a cube file by mp2 method.


Yi Yao,
1) Your input file is computing the electron density of the HF wavefunction, since MP2 is a perturbative method
that does not generate a new wavefunctions.
2) Could you please provide the results of your integration Python script?

Cheers, Edo

Just Got Here
4:26:46 PM PDT - Wed, Oct 16th 2013
My result is 41.8292212141 which should be 28 (11+17) in my opinion. Even it is the HF wavefunction, it should be 28. I double checked it by the bader analysis tool it give me exact the same answer.

So, is there some method to generate the charge density by some post-HF method like MP or CC?

My nwchem version is 6.1.1 on hopper at nersc.

Thank you very much,
Yi Yao

Forum Vet
10:12:31 PM PDT - Thu, Oct 17th 2013
Yi Yao

1) I have modified the dplot code so that the integrated density is printed now. In your NACal input,
the integrated value is close to the expected 28.
The patch can be applied to 6.3 (possibly to 6.1.1) and can be find at
http://nwchemgit.github.io/images/Dplot_intsum.patch.gz

2) There is a way to compute property from MP2, but it requires the calculation of MP2 gradients.
More details can be found in the documentation at
http://nwchemgit.github.io/index.php/Release62:MP2#One-electron_properties_and_natural_orbit...

Your input needs to be changes as followin
memory 1000 mb
start nacl

geometry
 na 0 0 0
 cl 2.0 0 0
end

basis spherical
 * library 6-31G*
end

mp2
freeze atomic
end

task mp2 gradient

dplot
 TITLE CHARGE_DENSITY
 vectors nacl.mp2nos
  LimitXYZ
-3. 5. 240
-3. 3. 180
-3. 3. 180
 spin total
 gaussian
 output naclmp2.cube
end

task dplot

Just Got Here
5:32:55 AM PDT - Fri, Oct 18th 2013
Thank you for your reply. It helps a lot!

And, I found the problem I got strange number of electrons.

It's the mesh I used to generate the charge density is not fine enough to capture the wiggling of wavefunction in some place especially for the inner shell electrons. When I change to plot the valence charge density and using a finer mesh, the problem solved.

Thank you for your help!

YY


Forum >> NWChem's corner >> Running NWChem