From NWChem
Viewed 3518 times, With a total of 4 Posts

Just Got Here
Threads 1
Posts 3


7: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 6311G**
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
 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])
 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()




Edoapra Forum:Admin, Forum:Mod, bureaucrat, sysop


Forum Vet
Threads 9
Posts 1431


1:03:39 PM PDT  Wed, Oct 16th 2013 

Quote:Yiy Oct 16th 7:42 amDear 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


3: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 postHF method like MP or CC?
My nwchem version is 6.1.1 on hopper at nersc.
Thank you very much,
Yi Yao




Edoapra Forum:Admin, Forum:Mod, bureaucrat, sysop


Forum Vet
Threads 9
Posts 1431


9: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.nwchemsw.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.nwchemsw.org/index.php/Release62:MP2#Oneelectron_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 631G*
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


4: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



AWC's:
2.5.10 MediaWiki  Stand Alone Forum Extension
Forum theme style by: AWC