edeprince3 / v2rdm_casscf

A variational 2-RDM-driven CASSCF plugin to Psi4

v2rdm_casscf

A variational 2-RDM-driven CASSCF plugin to Psi4

OVERVIEW

This plugin to Psi4[1] performs variational two-electron reduced-density-matrix (2-RDM)-driven complete active space self consistent field (CASSCF) computations. In principle, because variational 2-RDM (v2RDM) methods scale only polynomially with system size, vRDM-driven CASSCF computations can be performed using active spaces that are larger than can be used within conventional configuration-interaction-driven CASSCF methods. For more information regarding the performance of the method, see Refs. 2-4.

INSTALLATION

To run the Psi4 plugin v2rdm_casscf:

• Download Psi4 (1.2rc2 or later) from github.com: https://github.com/psi4/psi4, and follow the installation instructions given here: http://psicode.org/psi4manual/master/build_planning.html . Make sure to keep the name of the plugin directory v2rdm_casscf.

• Configure with CMake to generate a Makefile. Run psi4 --plugin-compile to get a CMake command. Modify it as needed with -D for compiler, libraries, and options.

• Note that, if you configured Psi4 with a fortran compiler, you shouldn't have to specify these things here. If the configure shows no errors, compile the plugin:

  make

• If the plugin compiles without any errors, you can run a few tests:

  cd tests

  make

• The test directories (tests/v2rdm1, etc.) contain input files that can help you get started using v2rdm-casscf.

• Or run tests through pytest

  pytest -rws -v tests/

  python -c "import v2rdm_casscf; v2rdm_casscf.test('quick')"

INPUT OPTIONS

N-representability conditions

• POSITIVITY (string):

  The positivity conditions enforced in the computation. Allowed values include DQG, DQ, DG, D, DQGT1, DQGT2, and DQGT1T2. The default value is DQG.

• CONSTRAIN_D3 (bool):

  Enforce the additional condition that D3 be possitive and correctly contract to D2? Default false.

• CONSTRAIN_SPIN (bool):

  Do constrain the expectation value of spin squared? Default true.

Convergence

• E_CONVERGENCE (double):

  The convergence in the primal/dual energy gap. Default 1e-4.

• R_CONVERGENCE (double):

  The convergence in the primal and dual errors. Default 1e-4.

• MAXITER (int):

  The maximum number of outer iterations. Default 10000.

Active space specification

• FROZEN_DOCC (array):

  An array containing the number of frozen doubly-occupied orbitals per irrep. These orbitals are not included in the active space, nor are they optimized during the v2RDM-CASSCF procedure. This option trumps Psi4's NUM_FROZEN_DOCC and FREEZE_CORE options.

• FROZEN_UOCC (array):

  An array containing the number of frozen unoccupied orbitals per irrep. These orbitals are not included in the active space, nor are they optimized during the v2RDM-CASSCF procedure. This option trumps Psi4's NUM_FROZEN_UOCC option.

• RESTRICTED_DOCC (array):

  An array giving the number of restricted doubly-occupied orbitals per irrep. These orbitals are not included in the active space, but they are optimized by the v2RDM-CASSCF procedure.

• RESTRICTED_UOCC (array):

  An array giving the number of restricted unoccupied orbitals per irrep. These orbitals are not included in the active space, but they are optimized by the v2RDM-CASSCF procedure.

• ACTIVE (array):

  An array giving the number of active orbitals (occupied plus unoccupied) per irrep. This option provides a more intuitive way of specifying the active space than the FROZEN_DOCC/RESTRICTED_DOCC/RESTRICTED_UOCC/FROZEN_UOCC keywords. The simplest specification of the active space would involve this keyword and at the RESTRICTED_DOCC keyword. This option trumps the RESTRICTED_UOCC option, which will be determined from the ACTIVE, RESTRICTED_DOCC, FROZEN_DOCC, and FROZEN_UOCC arrays.

Restarting jobs

• WRITE_CHECKPOINT_FILE (bool):

  Do save progress in a checkpoint file? Default false.

• CHECKPOINT_FREQUENCY (bool):

  Frequency of checkpoint file generation. The checkpoint file is updated every CHECKPOINT_FREQUENCY iterations. The default frequency will be ORBOPT_FREQUENCY.

• RESTART_FROM_CHECKPOINT_FILE (string):

  File containing previous primal/dual solutions and integrals.

Integrals and SCF type

• DF_BASIS_SCF (string):

  Auxiliary basis set for SCF density fitting computations. Defaults to a JKFIT basis.

• SCF_TYPE (string):

  What algorithm to use for the initial SCF computation. Default DF.

• CHOLESKY_TOLERANCE (double):

  Tolerance for Cholesky decomposition of the ERI tensor. Default 1e-4.

Orbital optimization

• ORBOPT_ONE_STEP (int):

  Flag to optimize orbitals using a quasi one-step type approach. Default 1.

• ORBOPT_GRADIENT_CONVERGENCE (double):

  Convergence in the orbital gradient norm. Default 1e-4.

• ORBOPT_ENERGY_CONVERGENCE (double):

  Convergence in the energy for orbital rotations. Default 1e-8.

• ORBOPT_EXACT_DIAGONAL_HESSIAN (int):

  Flag for using exact expresions for diagonal Hessian elements. Default 0.

• ORBOPT_FREQUENCY (int):

  Frequency of orbital optimization. Optimization occurs every ORBOPT_FREQUENCY iterations. Default 200.

• ORBOPT_ACTIVE_ACTIVE_ROTATIONS (bool):

  Do rotate active/active orbital pairs? Default false.

Additional files

• MOLDEN_WRITE (bool):

  Do write a MOLDEN output file containing the natural orbitals? If yes, the filename will end in .molden, and the prefix is determined by WRITER_FILE_LABEL (if set), or else by the name of the output file plus the name of the current molecule. Default false.

• WRITER_FILE_LABEL (string):

  Base filename for text files written by PSI, such as the MOLDEN output file, the Hessian file, the internal coordinate file, etc. Use the add_str_i function to make this string case sensitive.

• ORBOPT_WRITE (bool):

  Do write a ORBOPT output file? If so, the filename will end in .orbopt, and the prefix is determined by WRITER_FILE_LABEL (if set), or else by the name of the output file plus the name of the current molecule.

KNOWN ISSUES

• For large jobs, when running with multiple threads, sometimes a thread will hang and the job will stall.
• For large jobs, add "ulimit -s unlimited" to .bashrc to avoid segfault when calling the fortran orbital optimization routines.

REFERENCES

[1] J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. A. Evangelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke, M. L. Abrams, N. J. Russ, M. L. Leininger, C. L. Janssen, E. T. Seidl, W. D. Allen, H. F. Schaefer, R. A. King, E. F. Valeev, C. D. Sherrill, and T. D. Crawford, WIREs: Comp. Molec. Sci. 2, 556 (2012). "Psi4: an open-source ab initio electronic structure program"

[2] J. Fosso-Tande, D. R. Nascimento, and A. E. DePrince III, Mol. Phys. 114, 423-430 (2015). "Accuracy of two-particle N-representability conditions for describing different spin states and the singlet-triplet gap in the linear acene series." http://dx.doi.org/10.1080/00268976.2015.1078008

[3] J. Fosso-Tande, T.-S. Nguyen, G. Gidofalvi, and A. E. DePrince III, J. Chem. Theory Comput., 12, 2260-2271 (2016). "Large-scale variational two-electron reduced-density-matrix-driven complete active space self-consistent field methods." http://dx.doi.org/10.1021/acs.jctc.6b00190

[4] E. Maradzike, G. Gidofalvi, J. M. Turney, H. F. Schaefer III, and A. E. DePrince III, J. Chem. Theory Comput., 13, 4113-4122 (2017). "Analytic Energy Gradients for Variational Two-Electron Reduced-Density-Matrix-Driven Complete Active Space Self-Consistent Field Theory." http://dx.doi.org/10.1021/acs.jctc.7b00366

v2rdm_casscf's Issues

allow users to pick active space by reordering orbitals

Users should be able to specify total number of active occupied and active unoccupied orbitals. If they don't want the frontier orbitals, there should be a reorder array, like IORDER in gamess.

stop allocating memory on the stack in the fortran code

The orbital optimization code is allocating some huge stuff on the stack, making it necessary for users to set "ulimit -s unlimited" .

do not store core density

We're currently storing and passing the core-core part of the 2-RDM to the orbital optimization routines.

9-acene (dqg@cc-pvdz) crashes with seg. fault

tail of ouput before crash
    197   328 -1449.857504 -1869.967403  420.109899   1.000    2.50967    0.08929
    198   352 -1449.865533 -1867.154067  417.288534   1.000    2.51061    0.08912
    199   327 -1449.873647 -1864.346737  414.473090   1.000    2.51113    0.08895

    ==> Orbital Optimization <==

more test cases

We need more test cases. Model them off of v2rdm_casscf/tests/v2rdm1

write spin-free 1- and 2-rdm to output file

Add flag to write spin-free 1- and 2-RDM to output file.

additional tests

we need more/different tests. we need open-shell tests more than anything

Execution error ?

Hello,

When trying to compile the code with psi4public, the test example in the source dir, dump the following error,

An error has occurred python-side. Traceback (most recent call last):

File "", line 23, in

ImportError: ./v2rdm_casscf.so: undefined symbol: _ZN3psi6OEPropC2Ev

any idea, why is this happening ?

Thanks and best regards,
Jose Luis

print orbitals

Add flag to print orbitals.

A test project with v2rdm_casscf

Hi,

I am curious how v2RDM-CASSCF works and I did a test project inspired by the work of G. Chan.
http://aip.scitation.org/doi/pdf/10.1063/1.3275806

Basically, I want to calculate the isomerization energy of two isomers of Cu2O2. The active space is 28 electrons in 32 orbitals. The basis set is small. The molecule has D2h symmetry. Thus I think that it is feasible.
The reference result from DMRG is 25.6 kcal/mol.

This is an input file

sys.path.insert(0, '/temp0/tm/source')
import v2rdm_casscf
memory 120 Gb

molecule bis {
Cu 0.00  0.00 0.00
Cu 0.00  0.00 2.80
O  0.00  1.15 1.40
O  0.00 -1.15 1.40

2 1
units ang
}

basis mybas{
  assign Cu stuttgart
  assign O  roos-ano-dz
}

set globals{
  basis mybas
  reference rhf
  scf_type cd
  df_ints_io     save
  cholesky_tolerance 1.0e-4
#             Ag   B1g   B2g   B3g    Au   B1u   B2u   B3u 
  restricted_docc [ 4, 1, 0, 1, 0, 1, 2, 3]
  active          [ 6, 4, 4, 2, 2, 4, 6, 4]
       
  mcscf_maxiter 200
  nat_orbs true
  diis_max_vecs 20
  die_if_not_converged false
  print_mos true
}

scf_e, wfn = energy('b3lyp', return_wfn=True)
molden(wfn, 'b3lyp.molden')

set v2rdm_casscf {
  positivity                DQG
  r_convergence             1e-4
  e_convergence             1e-4
  maxiter                   10000
  constrain_d3              false
  write_checkpoint_file     true
  cholesky_tolerance        1e-4
  molden_write              true
  orbopt_write              true
}

cas_e, cas_wfn = energy('v2rdm-casscf', return_wfn=True)

Since maxiter is only 10000, the calculation crashed because v2RDM-CASSCF not converge.
The intriguing thing is that at the iteration 10000, the isomerization energy is only -0.8 kcal/mol. This error is huge and I don't know what causes this:

• It can be the poor initial active space, but I cannot find a way to manually choose it, by e.g. reading DFT or HF orbitals and rotate them. Is there any keyword to do this?
• In the end, there is no way I can check how the active space looks like.
• DQG+T2 should produce better results, but it is extremely slow and I gave up this calculation.

Thanks!

T1 and T2 are way too slow

The T1 and T2 conditions need to be rewritten in a more intelligent way.

rename jacobi routines and options

We aren't using the jacobi optimizer anymore. The options and functions should reflect that.

large jobs with multiple threads hang (sometimes...maybe)

sometimes the code hangs, but i haven't found a reproducible example yet.

make v2rdm_casscf compatible with new plugin framework

threaded orbital optimization

The orbital optimization should use the same number of threads as the rest of the code. There is no need for a separate input option.

LICENSE says GPLv3, but copyright boilerplate in source say GPLv2

Minor thing, but there's a discrepancy between the content of LICENSE and the (auto-generated) copyright headers in the source files: the former says GPLv3, while the latter say GPLv2+. Technically GPLv2 (or any later version) includes GPLv3, but it might make sense to clarify that.

flags to print 2-, 3-RDM

Other modules might eventually want to use the 2-RDM (or 3-RDM). We need input options and code to do this.