Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

CAS PES crashing #758

Open
susilehtola opened this issue Jul 6, 2017 · 5 comments

Comments

Projects
None yet
5 participants
@susilehtola
Copy link
Member

commented Jul 6, 2017

I'm trying to calculate a potential energy surface for H2O double dissociation in a (4e,4o) active space using the input

molecule h2o {
0 1
symmetry c1
units bohr
O
H 1 R
H 1 R 2 110.6
}

R0=1.84345

set basis cc-pVTZ
set scf_type pk
h2o.R=3*R0
scf_e, scf_wfn = energy('scf', return_wfn=True)

set frozen_docc [ 3 ]
set active [ 4 ]
set qc_module detci
set r_convergence 1e-8
set e_convergence 1e-10
set ci_maxiter 1000
set mcscf_algorithm ah

for ir in range(30, 7, -1):
    h2o.R=R0*ir/10
    print(h2o.R)
    if ir == 30:
       cas_e, cas_wfn = energy('casscf', ref_wfn=scf_wfn, return_wfn=True)
    else:
       cas_e, cas_wfn = energy('casscf', ref_wfn=cas_wfn, return_wfn=True)

Judging from output.dat, the run appears to crash in DiskJK for the third step.

Somewhat oddly, if I look at the occupation numbers for the second step, I see

   Active Space Natural occupation numbers:

         A   1620.561712         A   1145.170647         A   761.060310
         A   410.061564

whereas in the first step they're still sane

   Active Space Natural occupation numbers:

         A   1.119206         A   1.083963         A   0.916038
         A   0.880793

What's going on here?!

@fevangelista

This comment has been minimized.

Copy link
Contributor

commented Jul 7, 2017

This might be due to a problem we noticed when computing PES with CASSCF. Guess orbitals are assumed to be orthogonal, which is not true if the geometry is changed. We have developed some python scripts to manually orthogonalize the guess orbitals with respect to the current overlap matrix. @lcyyork may be able to help.

@lcyyork

This comment has been minimized.

Copy link

commented Jul 7, 2017

Here is a script can orthogonalize orbitals between different geometries. However it does not consider frozen orbitals. They can be simply added following the comments inside the script.

#!/usr/bin/python
import psi4

# some global variables
nirrep = None
nrdoccpi = None
nruoccpi = None
nactvpi = None
nmopi = None

"""
  This function makes (C1)^T S2 C1 orthogonal
  C1: converged CASSCF orbitals at geometry 1
  S2: SO overlap matrix at geometry 2
  return: orthogonal orbitals
"""
def ortho_orbs(wfn1, wfn2, semi = True):
    title = "\n  ==> Orthogonalize Orbitals Between Different Geometries <==\n"
    psi4.core.print_out(title)

    # make sure there is no frozen orbitals
    psi4.core.print_out("\n    Testing frozen orbitals ... ")
    global nirrep
    nirrep = wfn2.nirrep()
    nfdoccpi = psi4.core.Dimension.from_list(psi4.core.get_option("DETCI","FROZEN_DOCC"))
    nfuoccpi = psi4.core.Dimension.from_list(psi4.core.get_option("DETCI","FROZEN_UOCC"))
    nf = nfdoccpi.n() + nfuoccpi.n()
    if nf != 0:
        psi4.core.print_out("False")
        raise ValueError("I am too lazy to consider frozen orbitals.")
    else:
        psi4.core.print_out("Pass")

    # get C1 and S2
    C1 = wfn1.Ca()
    S2 = wfn2.S()

    # figure out irreps and orbital spaces
    global nmopi
    global nrdoccpi
    global nactvpi
    global nruoccpi
    nmopi = wfn2.nmopi()
    nrdoccpi = psi4.core.Dimension.from_list(psi4.core.get_option("DETCI","RESTRICTED_DOCC"))
    nactvpi = psi4.core.Dimension.from_list(psi4.core.get_option("DETCI","ACTIVE"))
    nruoccpi = psi4.core.Dimension(nirrep)
    for i in range(nirrep):
        nruoccpi[i] = nmopi[i] - nrdoccpi[i] - nactvpi[i]

    # create subspace orbitals: core, active, virtual
    psi4.core.print_out("\n    Preparing orbitals of subspaces ... ")
    Ccore = psi4.core.Matrix("C core", nmopi, nrdoccpi)
    Cactv = psi4.core.Matrix("C actv", nmopi, nactvpi)
    Cvirt = psi4.core.Matrix("C virt", nmopi, nruoccpi)

    # fill in data to orbitals of subspaces
    for h in range(nirrep):
        offset1 = nrdoccpi[h]
        offset2 = nactvpi[h] + offset1

        for i in range(nmopi[h]):
            # core
            for j in range(nrdoccpi[h]):
                Ccore.set(h, i, j, C1.get(h, i, j))

            # active
            for j in range(nactvpi[h]):
                Cactv.set(h, i, j, C1.get(h, i, j + offset1))

            # virtual
            for j in range(nruoccpi[h]):
                Cvirt.set(h, i, j, C1.get(h, i, j + offset2))
    psi4.core.print_out("Done")

    # orthogonalize core
    psi4.core.print_out("\n    Orthogonalizing orbitals of subspaces ... ")
    Ccore = ortho_subspace(Ccore, S2)
    
    # orthogonalize active
    Cactv = projectout(Cactv, Ccore, S2)
    Cactv = ortho_subspace(Cactv, S2)

    # orthogonalize virtual
    Cvirt = projectout(Cvirt, Ccore, S2)
    Cvirt = projectout(Cvirt, Cactv, S2)
    Cvirt = ortho_subspace(Cvirt, S2)
    psi4.core.print_out("Done")

    # fill in data to the new combined orbitals
    psi4.core.print_out("\n    Combining orbitals of subspaces ... ")
    Cnew = psi4.core.Matrix("new C", C1.rowdim(), C1.coldim())
    for h in range(nirrep):
        offset1 = nrdoccpi[h]
        offset2 = nactvpi[h] + offset1

        for i in range(nmopi[h]):
            # core
            for j in range(nrdoccpi[h]):
                Cnew.set(h, i, j, Ccore.get(h, i, j))

            # active
            for j in range(nactvpi[h]):
                Cnew.set(h, i, j + offset1, Cactv.get(h, i, j))

            # virtual
            for j in range(nruoccpi[h]):
                Cnew.set(h, i, j + offset2, Cvirt.get(h, i, j))
    psi4.core.print_out("Done")

    if semi:
        psi4.core.print_out("\n    Semicanonicalizing orbitals ... ")
        U = semicanonicalize(wfn2.Fa(), Cnew)
        Cnew = psi4.core.Matrix.doublet(Cnew, U, False, False)
        psi4.core.print_out("Done")

    psi4.core.print_out("\n\n")
    return Cnew

"""
  This function project CP out of C
  C: orbitals to be projected by P
  CP: orbitals of the projector
  S: SO overlap matrix
  return: Cp = (1 - P) C = C - CP (CP^T S C)
"""
def projectout(C, CP, S):
    M = psi4.core.Matrix.triplet(CP, S, C, True, False, False)
    P = psi4.core.Matrix.doublet(CP, M, False, False)
    Cp = C.clone()
    Cp.subtract(P)
    return Cp

"""
  This function orthogonalize C
  C: the orbitals of a subspace at geometry 1
  S: the SO overlap matrix at geometry 2
  return: C X where X is the canonical orthogonalizing transformation matrix
"""
def ortho_subspace(C, S):
    M = None
    try:
        M = psi4.core.Matrix.triplet(C, S, C, True, False, False)
    except ValueError as e:
        print "The dimensions are wrong. Cannot do Matrix::triplet."
    X = canonicalX(M)

    Cnew = psi4.core.Matrix.doublet(C, X, False, False)
    return Cnew

"""
  This function returns the canonical orthogonalizing transformation matrix
  S: overlap matrix
  return: X = U s^(-1/2)
"""
def canonicalX(S):
    rdim = S.rowdim()
    evals = psi4.core.Vector("evals", rdim)
    evecs = psi4.core.Matrix("evecs", rdim, rdim)

    S.diagonalize(evecs, evals, psi4.core.DiagonalizeOrder.Descending)
    shalf_inv = psi4.core.Matrix("s^(-1/2)", rdim, rdim)
    for h in range(nirrep):
        for i in range(rdim[h]):
            shalf_inv.set(h, i, i, evals.get(h, i) ** -0.5)

    X = psi4.core.Matrix.doublet(evecs, shalf_inv, False, False)
    return X

"""
  This function semicanonicalize the orbitals
  Fso: SO Fock matrix at geometry 2
  C: molecular orbitals that transforms Fso to Fmo
  return: unitary matrix that transforms orbitals to semicanonical orbitals
"""
def semicanonicalize(Fso, C):
    # transform SO Fock to MO Fock
    Fmo = psi4.core.Matrix.triplet(C, Fso, C, True, False, False)

    offsets = psi4.core.Dimension.from_list([0 * i for i in range(nirrep)])

    U = psi4.core.Matrix("U to semi", nmopi, nmopi)

    # diagonalize each blcok of Fmo
    for block in [nrdoccpi,nactvpi,nruoccpi]:
        F = psi4.core.Matrix("Fock",block,block)
        for h in range(nirrep):
            offset = offsets[h]
            for i in range(block[h]):
                for j in range(block[h]):
                    F.set(h, i, j, Fmo.get(h, i + offset, j + offset))
      
        evals = psi4.core.Vector("F Evals", block)
        evecs = psi4.core.Matrix("F Evecs", block, block)
        F.diagonalize(evecs, evals, psi4.core.DiagonalizeOrder.Ascending)
      
        for h in range(nirrep):
            offset = offsets[h]
            for i in range(block[h]):
                for j in range(block[h]):
                    U.set(h, i + offset, j + offset, evecs.get(h, i, j))
            offsets[h] += block[h] ### important ###

    return U

To use it, you can put the following to the input:

h2o.R = old
Ecas, wfn = energy('casscf', return_wfn=True)

h2o.R = new
Escf, wfnSCF = energy('scf', return_wfn=True)

# import the script <== Change HERE
sys.path.insert(0, '......')
from ...... import ortho_orbs

wfnSCF.Ca().copy(ortho_orbs(wfn,wfnSCF))
Ecas, wfn = energy('casscf', ref_wfn=wfnSCF, return_wfn=True)

I do not worry about efficiency usually, but you can probably make it faster using slicing when filling in the data to matrix.

@dgasmith

This comment has been minimized.

Copy link
Member

commented Jul 12, 2017

Anyone have a good idea where this function could live. Could be useful to have as the base part of Psi4.

@fevangelista

This comment has been minimized.

Copy link
Contributor

commented Jul 12, 2017

I think it should be called any time a guess set of orbitals is passed because there is no guarantee that MOs come from the same geometry and use the same overlap (metric) matrix. The precise spot, however, I do not know.

@loriab

This comment has been minimized.

Copy link
Member

commented Jul 12, 2017

Haven't vetted fully but should it be a py extension to Wfn? I've been thinking extensions should be better organized in driver.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.