[molpro-user] icpks convergence problem

[Tatiana Korona]

Hi,

Note that for E(20)ind,resp you don't need to do a finite field 
calculation, since you can calculate E(20)ind,resp as a 2nd order SCF 
property.

It is enough to save the effective SCF potential of monomer A, calculated 
as in the example script of Andreas, by adding in matrop a line:

load,Veffa,1211.2,square
save,Veffa,oper
}

and then to make for the monomer B:

{hf,thrcphf=1.e-10
polarizability,VEFFA
}

BTW, if the user-defined operator saved in this way and you still want to 
do the FF calculation, you add the field VEFFA simply by

field,VEFFA,0.001

No explicit loading/changing H0&H01 is necessary!

Best wishes,

Tatiana

Dr. Tatiana Korona http://tiger.chem.uw.edu.pl/staff/tania/index.html
Quantum Chemistry Laboratory
University of Warsaw
Pasteura 1, PL-02-093 Warsaw, POLAND


`The man who makes no mistakes does not usually make anything.'
                                        Edward John Phelps (1822-1900)

On Wed, 25 Apr 2007, Andreas Hesselmann wrote:

> Dear Glen,
>
> the default convergence tolerance for CPKS
> is set to 1d-6 for the density matrix.
> I have submit a patch (sapt_icpks) with which
> you can adjust this value (sapt_cpksthr) as
> well as the maximum number of iterations
> in the CPKS (which by default is set to 50)
> in the input.
>
> I would be sceptical about the result of a
> CPKS calculation if the convergence criteria
> of 1d-6 is not met within 50 iterations.
> You might check your result for the induction energy
> by comparing it with the finite field result.
>
> Here is an example for a finite field calculation
> of the induction energy for the water dimer:
>
>
> gthresh,energy=1d-11,orbital=1d-12
> geometry={nosym;noorient;angstroms;
> O1,, 0.0000000000, 0.0000000000, 0.0000000000
> H1,, 0.0000000000, 0.0000000000, 0.9621100000
> H2,, 0.0000000000, 0.9368337299,-0.2465760535
> O2,, 0.0000000000, 2.7788154372,-0.9708520107
> H3,,-0.7632231375, 2.9252149352,-1.5405904257
> H4,, 0.7632231375, 2.9252149352,-1.5405904257}
> basis=avdz
> int
>
> field=[0,0.001,-0.001]
>
> !---Monomer A-----------------
> dummy,o2,h3,h4
>
> {hf; maxit,20; save,2151.2}
>
> {matrop
> load,d,den,2151.2
> load,epot,epot
> load,h0,h0
> coul,J,d
> add,v,1,epot,2,J
> save,v,1211.2,square}
>
> !---Monomer B-----------------
> int
> dummy,o1,h1,h2
>
> {hf; start,atden; maxit,20; save,2152.2}
> e(1)=energy
>
> {matrop
> load,d,den,2152.2
> load,epot,epot
> load,h0,h0
> coul,J,d
> add,v,1,epot,2,J
> save,v,1212.2,square}
>
>
> delete,1
> !--->Induktion A->B
> {matrop
> load,v,square,1211.2
> load,h0,h0
> add,h01,h0,field(2),v
> save,h01,,h0}
>
> {hf; maxit,20; save,2153.2}
> e(2)=energy
>
> {matrop
> load,d,den,2153.2
> load,v,square,1211.2
> trace,v1,d,v
> set,v1,0.5*v1}
>
> delete,1
> {matrop
> load,v,square,1211.2
> load,h0,h0
> add,h01,h0,field(3),v
> save,h01,,h0}
>
> {hf; maxit,20; save,2154.2}
> e(3)=energy
>
> {matrop
> load,d,den,2154.2
> load,v,square,1211.2
> trace,v2,d,v
> set,v2,0.5*v2}
>
> show,v*
> table,e
> digits,12
> !values in mH
> e10_elst_A=-(e(2)-e(3))/(field(2)-field(3))*1000
> e20_ind_A= -0.5*(e(2)+e(3)-2*e(1))/((field(2)-field(1))*(field(3)-field(1)))*1000
> e20_ind_A=(v1-v2)/(field(2)-field(3))*1000
>
> text,Energies in mH
> table,e10_elst_A,e20_ind_A
> digits,12
>
> !--->Induktion B->A
> int
> dummy,o2,h3,h4
>
> {hf; start,atden; maxit,20; save,2155.2}
> e(1)=energy
>
> delete,1
> {matrop
> load,v,square,1212.2
> load,h0,h0
> add,h01,h0,field(2),v
> save,h01,,h0}
>
> {hf; maxit,20; save,2156.2}
> e(2)=energy
>
> {matrop
> load,d,den,2156.2
> load,v,square,1212.2
> trace,v1,d,v
> set,v1,0.5*v1}
>
> delete,1
> {matrop
> load,v,square,1212.2
> load,h0,h0
> add,h01,h0,field(3),v
> save,h01,,h0}
>
> {hf; maxit,20; save,2157.2}
> e(3)=energy
>
> {matrop
> load,d,den,2157.2
> load,v,square,1212.2
> trace,v2,d,v
> set,v2,0.5*v2}
>
>
> table,e
> digits,12
> !values in mH
> e10_elst_B=-(e(2)-e(3))/(field(2)-field(3))*1000
> e20_ind_B= -0.5*(e(2)+e(3)-2*e(1))/((field(2)-field(1))*(field(3)-field(1)))*1000
> e20_ind_B=(v1-v2)/(field(2)-field(3))*1000
> e20_ind=e20_ind_A+e20_ind_B
>
> text,Energies in mH
> table,e20_ind_A,e20_ind_B,e20_ind
> digits,8,8,8
>
>
> Of course the convergence thresholds for the SCF
> should be even higher in the finite-field case.
>
> Best regards,
> Andreas
>
>
>
>
>
>
>
>
> --------------------------------------------------
> Andreas Hesselmann
> Institut für Physikalische und Theoretische Chemie
> Universität Erlangen
> Egerlandstraße 3
> 91058 Erlangen / Germany
> Phone:  +49 9131/85-25021
> E-Mail: andreas.hesselmann at chemie.uni-erlangen.de
> -------------------------------------------------
>