How to add counterpoise (BSSE) correction for N₂–Rb interaction energy

[Joseph Van Forst]

Hi, 

                I'm  a beginner in MOLPRO and I'm trying to calculate the interaction energy of N2 and Rb. Here is my input file:

```

***,N2Rb

rhcl=[5.0,5.25,5.5,5.75,6.0,6.25,6.5,6.75,7.0,7.25,7.5,7.75,8.0,8.25]

do i=1,#rhcl
  r=rhcl(i)

geometry={
    N
    N,1,2.0775
    Rb,2,r,1,180.0
}

basis=aug-cc-pVQZ, rb=cc-pV5Z-PP


if(i.eq.1) then
    hf
  else
    hf
  end if

  mp2
  ccsd
  ccsd(t)
  eccsdt(i)=energy
end do

{table,rhcl,eccsdt
 plot,file='n2rb.plot'
 title, 'Results for N2Rb'
}

```

Is this the right approach to do this in MOLPRO?
How can I apply the counterpoise correction (Boys–Bernardi BSSE) to each point in the scan? Is there any direct way to do this? 

Thank you for your suggestion and help

[Peter Knowles]

Yes: https://www.molpro.net/manual/doku.php?id=intermolecular_interaction_energies

Peter

On 10 Aug 2025, at 12:22, Joseph Van Forst <rajg...@gmail.com> wrote:



External email to Cardiff University - Take care when replying/opening attachments or links.

Nid ebost mewnol o Brifysgol Caerdydd yw hwn - Cymerwch ofal wrth ateb/agor atodiadau neu ddolenni.

--
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, visit https://groups.google.com/d/msgid/molpro-user/84be84d5-f719-4ffe-b003-da1b6c3037f4n%40googlegroups.com.

[Joseph Van Forst]

Thank you for the response and I appreciate your help.  I have went through the documentation and and tried running this calculation with the INTERACT  procedure. However, I'm  now facing the error:

```

 ? Error
 ? SPIN not possible
 ? The problem occurs in putocc
 GLOBAL ERROR fehler on processor   0

```

Could you please tell me what am I doing wrong?

here is my simplified input file for a single distant point:

```

***,N2Rb BSSE corrected interaction energies
memory,500,m

rhcl=[8]

basis={
  default,aug-cc-pVQZ
  Rb=cc-pV5Z-PP
}

mycalc={
  hf
  mp2
  {ccsd, nocheck}
  {ccsd(t),nocheck}
}

do i=1,#rhcl
  r=rhcl(i)

  geometry={
    N, 0.0, 0.0, -1.03875
    N, 0.0, 0.0, 1.03875
    Rb, 0.0, 0.0, r
  }


  INTERACT,DO_CP=1,DO_NOCP=1,PROC=mycalc,UNIT=KCAL/MOL,SCALE=0.82

  eccsdt_cp(i)=E_INTERACT_CCSD(T)_CP
  eccsdt_nocp(i)=E_INTERACT_CCSD(T)_NOCP
end do

table,rhcl,eccsdt_cp,eccsdt_nocp
 plot,file='n2rb_cp_vs_nocp.plot'
 title,'CCSD(T) CP vs NOCP interaction energies for N2–Rb'
```

Thank you for your time and help. 

[andreas...@gmail.com]

Dear Joseph,

your whole complex has one unpaired electron, so you need to do the calculation with a doublet spin wave function. In INTERACT one sets the spin for each monomer with spin=[mon1,mon2,...] and the whole complex then is assumed to have the sum of the spins of the monomers.

So I think the correct input would be:

INTERACT,DO_CP=1,DO_NOCP=1,PROC=mycalc,UNIT=KCAL/MOL,SCALE=0.82,spin=[0,1]

with spin=0 for N2 (1st monomer) and spin=1 for Rb (2nd monomer).

However, I noticed that NTERACT treats each atom of N2 separately (so the whole complex as a trimer) the reason likely being that Molpro uses Angstrom over Bohr units (which I do not fully understand why this happens, normally he only switches to Angstrom if the geometry is being entered as xyz coordinates). To fix this, insert 'bohr' somewhere at the beginning, like:

rhcl=[8]
bohr                                  !use bohr units globally
basis={
  default,aug-cc-pVDZ
  Rb=cc-pVDZ-PP
}

The non-CP and CP interaction energies are collected in the variables DENOCP and DECP, respectively. In total Molpro saves them for 5 different methods in your example, see output (I used a smaller basis set for testing):

  Method               DE(NOCP)   DE(CP)
 HF                     1.561     1.607
 MP2                   60.805    61.002
 SCS-MP2               60.751    60.950
 CCSD                  61.209    61.395
 CCSD(T)               63.557    63.765

So if you want to just save the energies for CCSD(T) you can do this in the input by using:

 eccsdt_cp(i)=DECP(5)
 eccsdt_nocp(i)=DENOCP(5)

Best wishes,

Andreas

[hjwern...@gmail.com]

There are some issues here:

1.)  Andreas is right that BOHR needs to be given, since an xyz input is given and the default for this is ANG.

2.) It is not needed to specify MP2 and CCSD in the procedure. CCSD(T) does it all.

3.)  Rb has only 1 valence electron, and no correlation energy. For interaction energies with alkali atoms core-valence correlation effects are very important, and therefore core-correlation for Rb needs to be included. This is done by adding the directive core,mixed to CCSD(T). "Mixed" means that core correlation is included for elements on the left side of the periodic table, see manual. Of course, then a suitable basis set must be used (aug-cc-pwCVnZ).

4.) SCS-MP2 is not implemented for open-shell, and there is a bug in interact which causes this to be included in the second and following geometries of the loop over distances. This results in unreasonable SCS-MP2 interaction energies, which should be ignored. This can be avoided by using the variable option in interact (see attached input).

5.) The core-valence correlation effects are known to be very slowly convergent with basis set, and you should consider to use F12 or core-polarisation potentials to improve the results. Of course, double-zeta is too small in any case.  

6.) The interaction energies are stored in variable DECP and DENOCP.

Best regards

Joachim Werner

[Joseph Van Forst]

Dear Joachim Werner, 

               Thank you for your well structured and detailed answer. It is working fine. 

                                      As l learn more, I added  a "ghost" atom in the mid point of the Rb-N2 interaction vector  for adding "bond functions". When I add this "ghost atom" to the input file: the first step of the 'DO'  loop  ( r = 4.0 Bohr)is running fine. However in the second step, when (r = 9.0 Bohr), I get the following error:

"""

 MOLECULE CHANGED: INITIALIZING OCCUPATION AND DUMP RECORDS

 nelec=   8  ms2=   1

 ? Error
 ? SPIN not possible
 ? The problem occurs in putocc

 GLOBAL ERROR fehler on processor   0            
"""

This is confusing as the first stage is working just fine and the error is appearing only at the second stage. Could you take a look at my modified input file and let me know what could be causing this/

"""
rhcl=[8,9,10,11,12,13,14,15,16,18,20]  !distances in bohr
r= 8 ! to avoid zero initilization error

bohr      !use bohr units globally

geometry={
  N, 0.0, 0.0, -1.03875
  N, 0.0, 0.0, 1.03875

  o, 0, 0, r/2 ! Dummy atom
  Rb, 0.0, 0.0, r
}
dummy, o;
basis={
  N=aug-cc-pVDZ;
  Rb=aug-cc-pwCVDZ-PP;
  s,o,0.9,0.3,0.1;
  p,o,0.9,0.3,0.1;
  d,o,0.6,0.2;
  f,o,0.6,0.2;
  g,o,0.35;
}

mycalc={
  hf
  {ccsd(t);core,mixed}
  e=energy
  $method=ccsd(t)
}

do i=1,#rhcl
  r=rhcl(i)
  interact,do_cp=1,do_nocp=1,proc=mycalc,spin=[0,1],variable=e,scale=0.8
  eccsdt_cp(i)=DECP
  eccsdt_nocp(i)=DENOCP
end do

{table,rhcl,eccsdt_nocp,eccsdt_cp
plot,file='n2rb.plot'}
"""
I appreciate your kind help.

I'm attaching the full output: