Dplot output charge density, total number of electron incorrect

From NWChem

Viewed 2211 times, With a total of 4 Posts
Jump to: navigation, search
Forum >> NWChem's corner >> Running NWChem

Just Got Here
Threads 1
Posts 3
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
Threads 3
Posts 855
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
Threads 1
Posts 3
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
Threads 3
Posts 855
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://www.nwchem-sw.org/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://www.nwchem-sw.org/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
Threads 1
Posts 3
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



Who's here now Members 0 Guests 1 Bots/Crawler 0

AWC's: 2.5.10 MediaWiki - Stand Alone Forum Extension
Forum theme style by: AWC
Retrieved from "http://www.nwchem-sw.org/index.php/Special:AWCforum/st/id1013/"