Setting point charge field in molpro

[Bin Han]

Dear all,

I am looking for a way in molpro to add point charge field, and I have tried “LATTICE” option, but I find the result in molpro differs with the result in nwchem. I don't know whether the method is correct or there is any other method can achieve this. And in my nwchem code, I use the bq setting:

bq
     -0.20622000      -2.31432000        1.78993000 -1
end

and in molpro, I use:

 lattice,INFILE=input.molpro

and my input.molpro file:

# comment
1
-0.20622,-2.31432,1.78993,-1,1

but the results differ, in molpro :

 !UHF STATE 1.1 Energy               -710.995047465193
  UHF One-electron energy           -1216.928006086127
  UHF Two-electron energy             512.001251677005
  UHF Kinetic energy                  336.792835160185
  UHF Nuclear energy                   -6.068293056071
  UHF Virial quotient                  -2.111075335457

in nwchem:

        Total SCF energy =   -708.167180151266
      One electron energy =  -1231.647960876112
      Two electron energy =    529.549074219577
 Nuclear repulsion energy =     -6.068293494731

The results without point charge setting are same, in molpro:

 !UHF STATE 1.1 Energy               -707.816596217750
  UHF One-electron energy           -1237.414547945169
  UHF Two-electron energy             529.597951727420
  UHF Kinetic energy                  352.992689577910
  UHF Nuclear energy                    0.000000000000
  UHF Virial quotient                  -2.005187691179

in nwchem:

         Total SCF energy =   -707.816514794644
      One electron energy =  -1237.413776532485
      Two electron energy =    529.597261737841
 Nuclear repulsion energy =      0.000000000000
 

In my opinion, I think the "LATTICE" setting is just like putting a point charge in the position which is the same function using "bq" in nwchem, but I can't understand the different results. So hope you can give me some suggestions.

Best regard,

Bin

[Peterson, Kirk]

Dear Bin,

 

I am not an expert at all in applying lattice charges, but I do note a couple of things:

 

1) your UHF in Molpro did not converge, so it's not a valid comparison

2) the nuclear repulsion energy is the same between NWChem and Molpro, so the lattice charges do seem to be identical in the two programs

 

Have you tried comparing UHF energies for this input without the lattice charges?  I note that NWChem complains a bit about the local part of the ECP on Eu.  If you could post a simpler input that exhibits problems it would also be helpful to debug.

 

regards,

 

-Kirk

--
You received this message because you are subscribed to the Google Groups "molpro-user" group.
To unsubscribe from this group and stop receiving emails from it, send an email to molpro-user...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/af102a06-7379-4fc8-a604-988e70ab5d47n%40googlegroups.com.

[Bin Han]

Dear Kirk,

Thanks for your reply. I add a shift in UHF calculation and I get the converged result:

 !UHF STATE 1.1 Energy               -711.336479663842
  UHF One-electron energy           -1233.376651724358
  UHF Two-electron energy             528.108465116586
  UHF Kinetic energy                  352.754911802277
  UHF Nuclear energy                   -6.068293056071
  UHF Virial quotient                  -2.016517575984

The result still differs from the result in nwchem. If we check the results carefully, we can find that in nwchem the two-electron energy with point charge changes very small comparing with the result without point charge, but in molpro, we can find that the two-electron and one-electron energy all change very large with the result without point charge. So I guess there may be some difference dealing with these two energies in the two programs. In nwchem, I find an example in the manual:

and in molpro, the introduction about "LATTICE":

And actually I have compared UHF energies for this input without the lattice charges, and the results are identical:

The results without point charge setting are same, in molpro:

 !UHF STATE 1.1 Energy               -707.816596217750
  UHF One-electron energy           -1237.414547945169
  UHF Two-electron energy             529.597951727420
  UHF Kinetic energy                  352.992689577910
  UHF Nuclear energy                    0.000000000000
  UHF Virial quotient                  -2.005187691179

in nwchem:

         Total SCF energy =   -707.816514794644
      One electron energy =  -1237.413776532485
      Two electron energy =    529.597261737841
 Nuclear repulsion energy =      0.000000000000

Peterson, Kirk <kipe...@wsu.edu> 于2022年3月12日周六 23:59写道：

[Bin Han]

Dear Kirk,

Thanks for your reply. I add a shift in UHF calculation and I get the converged result:

 !UHF STATE 1.1 Energy               -711.336479663842
  UHF One-electron energy           -1233.376651724358
  UHF Two-electron energy             528.108465116586
  UHF Kinetic energy                  352.754911802277

  UHF Nuclear energy                   -6.068293056071
  UHF Virial quotient                  -2.016517575984

The result still differs from the result in nwchem. If we check the results carefully, we can find that in nwchem the two-electron energy with point charge changes very small comparing with the result without point charge, but in molpro, we can find that the two-electron and one-electron energy all change very large with the result without point charge. So I guess there may be some difference dealing with these two energies in the two programs. In nwchem, I find an example in the manual:

and in molpro, the introduction about "LATTICE":

And actually I have compared UHF energies for this input without the lattice charges, and the results are identical:

The results without point charge setting are same, in molpro:

 !UHF STATE 1.1 Energy               -707.816596217750
  UHF One-electron energy           -1237.414547945169
  UHF Two-electron energy             529.597951727420
  UHF Kinetic energy                  352.992689577910
  UHF Nuclear energy                    0.000000000000
  UHF Virial quotient                  -2.005187691179

in nwchem:

         Total SCF energy =   -707.816514794644
      One electron energy =  -1237.413776532485
      Two electron energy =    529.597261737841
 Nuclear repulsion energy =      0.000000000000

Best regards,

Bin

[Peterson, Kirk]

Can you try a simpler example without an ECP center?  Perhaps just Ar atom?

 

regards,

 

-Kirk

 

From: Bin Han <ice86...@gmail.com>
Date: Saturday, March 12, 2022 at 6:39 PM
To: Peterson, Kirk <kipe...@wsu.edu>
Cc: molpro-user <molpr...@googlegroups.com>
Subject: Re: [molpro-user] Setting point charge field in molpro

Dear Kirk,

 

Thanks for your reply. I add a shift in UHF calculation and I get the converged result:

 

 !UHF STATE 1.1 Energy               -711.336479663842
  UHF One-electron energy           -1233.376651724358
  UHF Two-electron energy             528.108465116586
  UHF Kinetic energy                  352.754911802277
  UHF Nuclear energy                   -6.068293056071
  UHF Virial quotient                  -2.016517575984

 

The result still differs from the result in nwchem. If we check the results carefully, we can find that in nwchem the two-electron energy with point charge changes very small comparing with the result without point charge, but in molpro, we can find that the two-electron and one-electron energy all change very large with the result without point charge. So I guess there may be some difference dealing with these two energies in the two programs. In nwchem, I find an example in the manual:

 

and in molpro, the introduction about "LATTICE":

 

[Bin Han]

Deal all,

I have tried a simple case only including Ar atom at (0,0,0) with a negative point charge at (-0.2,0,0). The result in molpro still differ from the result in nwchem:

in molpro:

 !UHF STATE 1.1 Energy               -545.377521599308
  UHF One-electron energy            -690.876943101977
  UHF Two-electron energy             193.125370483938
  UHF Kinetic energy                  517.176792426756
  UHF Nuclear energy                  -47.625948981270
  UHF Virial quotient                  -1.054528218562
in nwchem:

         Total SCF energy =   -544.560394739786
      One electron energy =   -689.415450345573
      Two electron energy =    192.481008030378
 Nuclear repulsion energy =    -47.625952424590

It's quite confusing, so I guess the two programs may use different models dealing with the point charge? 

Best regards,

Bin

在2022年3月13日星期日 UTC+8 01:18:05<kipe...@wsu.edu> 写道：

[Bin Han]

I have tried the q-chem , the result is same as nwchem:

 Hartree-Fock
 A restricted SCF calculation will be
 performed using DIIS
 SCF converges when DIIS error is below 1.0e-05
 ---------------------------------------
  Cycle       Energy         DIIS error
 ---------------------------------------
    1    -532.9717191935      4.14e-01
    2    -543.3488926565      2.43e-01
    3    -544.5421044486      2.72e-02
    4    -544.5603155083      1.60e-03
    5    -544.5603876983      3.29e-04
    6    -544.5603921951      5.59e-05
    7    -544.5603924886      1.40e-05
    8    -544.5603924811      5.85e-07  Convergence criterion met
 ---------------------------------------
 SCF time:   CPU 0.10s  wall 0.00s
 SCF   energy in the final basis set =     -544.5603924811
 Total energy in the final basis set =     -544.5603924811

So why does molpro give a wrong result?

Bin Han <ice86...@gmail.com> 于2022年3月13日周日 17:51写道：

To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/5384a6a3-d36b-4a4a-bcee-aefec3857447n%40googlegroups.com.

[Peterson, Kirk]

I'm hoping someone who is more familiar with this part of the code can address this.

To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/CAEbcpXSrjYOzn%3D1Xq3-D-zCV7fAZ%3DrHMi3i2uLvpufsQ0daO_w%40mail.gmail.com.

[Klaus Doll]

There seems to be a bug when using the keyword "cartesian" together with point charges.
For the time being, do not use "cartesian" with point charges.

Best wishes
Klaus

To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/MWHPR01MB3231B2665B53556A4FBD0AF5D60E9%40MWHPR01MB3231.prod.exchangelabs.com.

[Hans-Joachim Werner]

Dear Bin,
I am sorry for the trouble. It turns out that the default lattice method in Molpro does not support cartesian basis functions. Please use

lattice,INFILE=input.molpro,method=argos !(alternatively method=xfdint or method=potint also work, but may be a bit slower)

which gives for your example -544.56039254. The default method can be an order of magnitude faster to compute the lattice integrals, but unfortunately it does not work with cartesians. Sorry that we did not notice this earlier, since we hardly ever use cartesians.

In general, for density fitting and explicitly correlated methods cartesian basis functions are not recommended. The integrals may be MUCH slower (if they work at all).

We will change the program to catch the cartesian case properly for lattice integrals.

Best wishes
Joachim

---
Prof. Dr. Hans-Joachim Werner
Institut für Theoretische Chemie
Universität Stuttgart
Pfaffenwaldring 55
70569 Stuttgart, Germany
e-mail: wer...@theochem.uni-stuttgart.de

> To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/CAEbcpXSrjYOzn%3D1Xq3-D-zCV7fAZ%3DrHMi3i2uLvpufsQ0daO_w%40mail.gmail.com.

[Bin Han]

Dear Joachim,

Thanks for your reply, this new method setting finally solve my problem. However,  the "cartesian" setting is very important for my calculation, so I may can not remove the setting. By the way, actually, I need to do the CASSCF and CASPT2 calculation next, so I'd like to consult you that whether the "cartesian" setting will influence the result of CASSCF and CASPT2 calculation? Since programs that can do CASPT2 calculation are not so much, I may can not use another program to do a benchmark or reference calculation, so I am not sure about the result of CASSCF and CASPT2 calculation with "cartesian" setting. As a result, it is appreciated about any comments for my calculation with "cartesian" setting.

Best regards,

Bin