Align the phase of wavefunctions at different geometry.

[Xuezhi Bian]

Dear Molpro community, 

I am trying to use MRCI method (or any other excited state methods) to produce a multistate potential energy surface and meanwhile I want to get matrix elements between different states (e.g., spin-orbit coupling or dipole moment). However, the surface of such matrix elements is never smooth since there is sign (gauge) of wavefunction is arbitrary. 

I am wondering if there is a way to align the phase between wavefunctions at nearby geometries or if there is a way to calculate the overlap between them so that I can align them myself. 

Thanks very much for your time,

Xuezhi 

[Peterson, Kirk]

Dear Xuezhi,

 

yes, the only way is to calculate the overlap between the states and manually keep track of the sign that way:

 

ci;trans,4100.2,4200.2,ov

 

where 4100.2 and 4200.2 contain MRCI wave functions.

 

regards,

 

-Kirk

 

From: <molpr...@googlegroups.com> on behalf of Xuezhi Bian <xzb...@sas.upenn.edu>
Date: Tuesday, January 24, 2023 at 12:27 AM
To: molpro-user <molpr...@googlegroups.com>
Subject: [molpro-user] Align the phase of wavefunctions at different geometry.

 

[EXTERNAL EMAIL]

--
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/e9ea8318-d2fb-4cac-bd66-1cceaf2400dcn%40googlegroups.com.

[Alexander MITRUSHCHENKOV]

You should rather use bi-orthogonal version if geometries are different, otherwise ci;trans does not take into account the orbital change. So please use:

ci;trans,4100.2,4200.2,biorth

Alexander

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

Dr. Alexander Mitrushchenkov, IGR
Laboratoire de Modélisation et Simulation Multi Echelle
UMR 8208 CNRS
Univ Gustave Eiffel
5 Bd Descartes
77454 Marne la Vallée, Cedex 2, France

Phone:    +33(0)160957316
Fax:      +33(0)160957320
e-mail:   Alexander.Mitrushchenkov AT univ-eiffel.fr

To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/E7D1EE6E-F192-4D92-B695-60B28E5F67E1%40wsu.edu.

[Peterson, Kirk]

Thanks Sasha.  Last time I did this (a long time ago) the biorth option wasn't available.

 

-Kirk

[Xuezhi Bian]

Dear Kirk and Alexander,

Thanks so much for your answers. I tried to follow your instructions, but I am confused by the results I got. 

The ci overlap computed from biorth method is always 0 in my calculations. There might be a bug in my code, so I attached the input and output files. 

Actually, I noticed that when I tried to save the ci record from last step to a file using "data, rename". The record is appended to the file instead of overwritten the record with the same name. 

I am not sure how molpro choose from these records in the same file with the same name.  This might be the bug and there must be smarter ways to do it.  I really appreciate it if anyone can shed some light on this.

Thanks again,

Xuezhi

[Peterson, Kirk]

Dear Xuezhi,

 

I'm not sure I'm following the logic in your input very well.  Why don't you just increment the record for the CI wavefunctions by your do loop counter?

 

e.g.:

 

{ci,CIREC=5100.2,NOCHECK,thrden=1e-5,thrvar=1e-4;print,ci=0;maxiter,200,2000;wf,elec=18,spin=0,sym=1;save,3010.2+i;noexc}

{ci,CIREC=5200.2,NOCHECK,thrden=1e-5,thrvar=1e-4;print,ci=0;maxiter,200,2000;wf,elec=18,spin=2,sym=1;save,3020.2+i;noexc}

{ci,CIREC=5300.2,NOCHECK,thrden=1e-5,thrvar=1e-4;print,ci=0;maxiter,200,2000;wf,elec=18,spin=2,sym=2;save,3030.2+i;noexc}

 if (i.gt.1) then

{ci; trans,3010.2+i-1,3010.2+i,biorth}

{ci; trans,3020.2+i-1,3020.2+i,biorth}

{ci; trans,3030.2+i-1,3030.3+i,biorth}

endif

To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/4ed5c804-555d-4ee2-85de-e257fd5e9d12n%40googlegroups.com.

[Xuezhi Bian]

Thanks Kirk,

This is really helpful. I didn't realize I could increment the record number before. 

Now it makes much more sense, and I get some meaningful numbers (usually close to 1 for nearby geometry) by calculating only ci coefficient overlaps (ci, trans, ov). 

But calculating biorthogonalization overlap still gives me 0.  I don't know how this method works, but I never save or give it the orbitals for corresponding geometries as input. 

Did I miss something?

Thanks again,

Xuezhi

[Peterson, Kirk]

Dear Xuezhi,

 

glad it is more or less working now.  For just keeping track of the signs, doing it without biorthogonalization should work fine I think.  As to why that is always giving zero, I will defer to Alexander since that is his code.

 

regards,

 

-Kirk

 

From: Xuezhi Bian <xzb...@sas.upenn.edu>

Date: Tuesday, January 24, 2023 at 1:59 PM
To: "Peterson, Kirk" <kipe...@wsu.edu>

[Alexander MITRUSHCHENKOV]

Dear Xuezhi,

As I said, you should normally never use ci;trans for different geometries without bi-orthogonalization as you have no guarantee that the orbitals preserve the phase. The only exception is if you check that orbitals are similar, or e.g. use diabatic orbitals which more or less ensure this.

Bi-orthogonalization won’t work if geometries are too far away, but in this case overlap is nearly zero anyway. So use smaller steps. Also, increasing basis set may be useful as diffuse functions overlap more.

In your input I’ve also find many unnecessary things. You note that now MULTI can save directly CI vectors into CI records, so no subsequent MRCI is needed; for different state symmetries the record number is increased by 100.

The following input works fine for me.

Best wishes,

Alexander

=====

symmetry,z

bohr

ir1=1.1

theta=89

i=0

do ir2=0.1,10.0,0.2

i=i+1

r1(i)=ir1

r2(i)=ir2

angle(i)=theta

xh1=-1*r1(i)/2

xh2=r1(i)/2

xs=r2(i)*cos(theta)

ys=r2(i)*sin(theta)

geometry={

3

h,xh1,0.0,0.0

h,xh2,0.0,0.0

s,xs,ys,0.0

}

if (i.gt.1) then

data,rename,3010.2,3011.2

data,rename,3110.2,3111.2

data,rename,3210.2,3211.2

endif

basis=vdz

{multi, cisave=3010.2

maxiter,50; occ, 13,2; closed, 4,1;

wf,elec=18,spin=0,sym=1;

state,1;

wf,elec=18,spin=2,sym=1;

state,1;

wf,elec=18,spin=2,sym=2;

state,1;

}

emrci1(i)=energy(1)

emrci2(i)=energy(2)

emrci3(i)=energy(3)

if (i.gt.1) then

{ci; trans,3010.2,3011.2,biorth}

{ci; trans,3110.2,3111.2,biorth}

{ci; trans,3210.2,3211.2,biorth}

endif

lsint

{ci; tranls, 3010.2,3110.2,0,0}

lz(i)=trlsz

{ci; tranls, 3010.2,3210.2,0,2}

lx(i)=trlsx

ly(i)=trlsy

enddo

{table,r1,r2,angle,emrci1,emrci2,emrci3,lx,ly,lz

save,H2S.dat

}

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

Dr. Alexander Mitrushchenkov, IGR
Laboratoire de Modélisation et Simulation Multi Echelle
UMR 8208 CNRS
Univ Gustave Eiffel
5 Bd Descartes
77454 Marne la Vallée, Cedex 2, France

Phone:    +33(0)160957316
Fax:      +33(0)160957320
e-mail:   Alexander.Mitrushchenkov AT univ-eiffel.fr

To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/022E37BD-20EE-429A-865D-B8D16EEABD80%40wsu.edu.

[Xuezhi Bian]

Dear Alexander,

Thank you very much for your help. I understand this much better now. 

I did some extra unnecessary things in my script because I naively thought the format of ci-vector record from multi is different from the format of record that mrci part need to be given as an input.  Thanks for clearing things up. 

I have one last little question. Is the overlap calculated from {ci, trans, biorth}  (i.e. <1.1||1.1>) assigned to any system recognized variable so that I can manipulate them in the same script?  

Thanks again,

Xuezhi

[Alexander MITRUSHCHENKOV]

Sure, this is TROV variable, rectangular storage with ket states running faster I guess.

Best wishes,

Alexander

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

Dr. Alexander Mitrushchenkov, IGR
Laboratoire de Modélisation et Simulation Multi Echelle
UMR 8208 CNRS
Univ Gustave Eiffel
5 Bd Descartes
77454 Marne la Vallée, Cedex 2, France

Phone:    +33(0)160957316
Fax:      +33(0)160957320
e-mail:   Alexander.Mitrushchenkov AT univ-eiffel.fr

To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/3a9c410e-e30e-4f90-b7b5-9fb499812b70n%40googlegroups.com.