[molpro-user] Trouble computing nonadiabatic couplings with DDR

[Robert Send]

Dear Molpro-users,

I am trying to compute nonadiabatic coupling vectors with the input 
script below, using the DDR method. However, depending on whether or 
not I compute the couplings at bond distance r=0.5 before the rest of 
the distances (1.0,2.0,3.0,4.0), I get different results. In other 
words,
r=[0.5,1.0,2.0,3.0,4.0]
gives different couplings from
r=[1.0,2.0,3.0,4.0].

I tried to stick as close as possible to the LiF input script from the 
molpro online manual.

Thanks a lot for your help,

Robert

Here comes the input:

memory,160000000
basis,h=aug-cc-pvdz     r=[1.0,2.0,3.0,4.0]
dr=0.001               geometry={
x,y
h
h,h,rbh
}

hf

do i=1,#r             rbh=r(i)          hf

    {casscf
    closed
    occ,10,4,4,0
    wf,2,1,0
    state,2,1,2
    orbital,2140.2
    }

   {ci;
   wf,2,1,0
   state,2,1,2   noexc            save,6000.2;    dm,8000.2}

   rbh=r(i)+dr

    {casscf
    closed
    occ,10,4,4,0
    wf,2,1,0
    state,2,1,2
    start,2140.2
    orbital,2141.2
    }

   {ci;
    wf,2,1,0
    state,2,1,2
    noexc;
    save,6001.2}

   {ci;trans,6000.2,6001.2;       dm,8100.2}

    rbh=r(i)-dr

   {casscf
    closed
    occ,10,4,4,0
    wf,2,1,0
    state,2,1,2
    start,2140.2
    orbital,2142.2
    }

   {ci;
    wf,2,1,0
    state,2,1,2
    noexc;save,6002.2}

   {ci;trans,6000.2,6002.2;         dm,8200.2}

   {ddr,dr,2140.2,2141.2,8100.2}    nacme1p(i)=nacme                 
{ddr,-dr,2140.2,2142.2,8200.2}   nacme1m(i)=nacme

   {ddr,2*dr                         orbital,2140.2,2141.2,2142.2;    
density,8000.2,8100.2,8200.2}     nacme2(i)=nacme

end do

nacmeav=(nacme1p+nacme1m)*0.5      
table,r,nacme1p,nacme1m,nacmeav,nacme2  title,Non-adiabatic couplings 
for H2

---------------------------------------

Robert Send
Theoretische Chemie
Institut für Physikalische Chemie
Universität Karlsruhe

+49-721-608-8344