Write input file for caspt2 optimization

huangtf · 2017-04-25 09:10:46

Dear team,

I am new to Molcas. Recently, I am reading the manual in the Molcas 8.0/doc/manual folder. I will try to run the calculation of CASPT2/ANO-RCC-VTZP (12,12) optimization for singlet NUN on the ground state that is an actinide molecule, and I have wrote an input file as shown below. However, I'm not sure whether it will be performed smoothly by serial Molcas 8.0 on the CentOS_6 system. Because I think my input file is incorrect.

Please help me out with your invaluable suggestions. Thank you!

Best regards,
huangtf

&gateway 
 title=NUN caspt2 optimization
 coord=NUN.xyz
 basis=ano-rcc-vtzp
 group=nosym

>>> export molcas_maxiter=500
>>> do while

 &seward; cholesky
 &rasscf; nactel=12 0 0; inactive=6
  ras2=12; lumorb
 &caspt2; frozen=41 
 &slapaf 

>>> enddo

&mckinley

Last edited by huangtf (2017-04-25 09:14:07)

Ignacio · 2017-04-25 09:56:45

I don't see anything wrong (provided the wavefunction calculation is correct), but note that all CASPT2 optimizations will use numerical gradients, so they may be slower than desired.

huangtf · 2017-04-25 16:09:18

Dear Ignacio,

Thank you very much for your reply!
What should I do to ensure that the wavefunction calculation is correct? Maybe this question is very childish, beg your pardon, please!

Thanks in advance!

Sincerely yours,
huangtf

Last edited by huangtf (2017-04-26 04:12:41)

Ignacio · 2017-04-26 08:24:19

Basically, check that the active after the CASSCF calculation is the one you wanted.

By the way, that "frozen=41" in CASPT2 looks suspicious. With 6 inactive orbitals and 12 active, freezing 41 orbitals looks like you are freezing your whole system and CASPT2 is actually doing nothing. More often than not you can just omit the "frozen" keyword.

huangtf · 2017-04-26 13:15:58

Dear Ignacio,

Thank you for your kind and helpful suggestion, I will delete the "frozen=41" in my input file. I still have three questions to consult you, as listed below.

1. Do you think some other modules, keywords or EMIL commands should be added into my input file to improve it?

2. At the beginning of the last paragraph of page123 in the manual, “The CASPT2 method can be used in any case where a valid reference function can be obtained with the CASSCF method.”, here, “a valid reference function can be obtained” means that the wavefunction calculation is correct, and in other word, the active after the CASSCF calculation is the one I wanted. Do you consider what I understand is right?

3. My last question is that, it seems impossible to check the active after the CASSCF calculation in time, because I feel the CASPT2 module will be performed immediately by Molcas after the CASSCF calculation. Do you think so?

Thanks a lot!

Best regards,
huangtf

Ignacio · 2017-04-26 17:20:15

huangtf wrote:

1. Do you think some other modules, keywords or EMIL commands should be added into my input file to improve it?

I would use the "RICD" keyword in &GATEWAY and remove "Cholesky" in &SEWARD, but that's mostly a matter of personal preference in this case.

2. At the beginning of the last paragraph of page123 in the manual, “The CASPT2 method can be used in any case where a valid reference function can be obtained with the CASSCF method.”, here, “a valid reference function can be obtained” means that the wavefunction calculation is correct, and in other word, the active after the CASSCF calculation is the one I wanted. Do you consider what I understand is right?

I understand the the CASSCF calculation should be converged. Even if the active space is "wrong" you can apply the CASPT2 method and get some energies, but they may not be the ones you want.

3. My last question is that, it seems impossible to check the active after the CASSCF calculation in time, because I feel the CASPT2 module will be performed immediately by Molcas after the CASSCF calculation. Do you think so?

You can do a single-point CASSCF calculation first:

&gateway 
 title=NUN single-point casscf
 coord=NUN.xyz
 basis=ano-rcc-vtzp
 group=nosym

 &seward; cholesky
 &rasscf; nactel=12 0 0; inactive=6
  ras2=12; lumorb

and check the result before starting the CASPT2 optimization. Or you can run the CASPT2 optimization and open the $Project.rasscf.molden file while it is running (but after it's finished &RASSCF).

huangtf · 2017-04-28 15:30:10

Dear Ignacio,

Thanks a lot for your reply last time! I have three questions to consult you.

1. By reading page 60 in the manual, I learn that the keyword “inactive” in &RASSCF specifies the number of doubly occupied orbitals (in each symmetry) that will not be included in the electron excitations and thus remain doubly occupied throughout the calculation, and these orbitals will be optimized. Do you agree with me?

2. According to “Another special case occurs if both CHARGE and NACTEL are given in the input and there is no symmetry, then the default value of INACTIVE will be automatically determined.” that was presented on line 7 from the bottom of page 294 in the manual or “Default is 0 in all symmetries, but if there is no symmetry (C1) and both CHARGE and NACTEL are given, the number of inactive orbitals will be calculated automatically.” seen in the paragraph describing the meaning of the keyword “inactive” in &RASSCF (page 300), for NUN(I) cation, supposing the active space is 11 electrons in 12 orbitals, no symmetry imposed, “frozen=30” in &RASSCF and without “inactive”, and then the number of inactive orbitals calculated automatically will be 17. Do you think so?

3. Do you think the inactive=“6” in &RASSCF of my input file (12 electrons in 12 orbitals) is reasonable or correct, and some other values of “inactive” are ok or better?

Thank you in advance!

Sincerely yours,
huangtf

Ignacio · 2017-04-28 16:00:34

huangtf wrote:

1. By reading page 60 in the manual, I learn that the keyword “inactive” in &RASSCF specifies the number of doubly occupied orbitals (in each symmetry) that will not be included in the electron excitations and thus remain doubly occupied throughout the calculation, and these orbitals will be optimized. Do you agree with me?

Yes, inactive orbitals are optimized. Only frozen (and deleted) orbitals are not optimized.

2. According to “Another special case occurs if both CHARGE and NACTEL are given in the input and there is no symmetry, then the default value of INACTIVE will be automatically determined.” that was presented on line 7 from the bottom of page 294 in the manual or “Default is 0 in all symmetries, but if there is no symmetry (C1) and both CHARGE and NACTEL are given, the number of inactive orbitals will be calculated automatically.” seen in the paragraph describing the meaning of the keyword “inactive” in &RASSCF (page 300), for NUN(I) cation, supposing the active space is 11 electrons in 12 orbitals, no symmetry imposed, “frozen=30” in &RASSCF and without “inactive”, and then the number of inactive orbitals calculated automatically will be 17. Do you think so?

NUN+ has 105 electrons (if it's nitrogen-uranium-nitrogen). 11 active electrons leaves 94 for 47 doubly-occupied orbitals, if 30 of these are frozen you should have 17 inactive orbitals. I'm not 100% sure the program takes into account the frozen orbitals, try and report back if it doesn't, because it should.

3. Do you think the inactive=“6” in &RASSCF of my input file (12 electrons in 12 orbitals) is reasonable or correct, and some other values of “inactive” are ok or better?

In your original input you have 12 active electrons and 12 (in 6 orbitals) inactive. That's a total of 24 electrons, so you would have a total charge of +82. I guess what you originally wanted was to put the "frozen=41" option in &RASSCF, not in &CASPT2.

As for the right number of inactive orbitals, it's given by the total number of electrons and the number of active electrons (which depends on the chemical properties and system you want to study, typically part of the valence electrons) and the number of frozen orbitals (typically some inner-core orbitals, if any)

huangtf · 2017-05-02 08:46:13

Dear Ignacio,

Thank you very much for your detailed reply last time! I have three questions to consult you.

1. “The default is to freeze the max of those that were frozen in the RASSCF calculation and the deep core orbitals.” in the paragraph describing the meaning of the keyword “frozen” in &CASPT2 (page 128 in the manual) implies that, without the “frozen” keyword in &CASPT2, the default value of “frozen” in &CASPT2 will be equal to the set value of “frozen” in &RASSCF when the latter is equal to or more than the number of the deep core orbitals, and will be more than the latter when the latter is less than the number of the deep core orbitals. To keep that the default value of “frozen” in &CASPT2 is equal to the set value of “frozen” in &RASSCF and to avoid that freezing too many orbitals looks like I am freezing my whole system and thus &CASPT2 is actually doing nothing, I should set the value of “frozen” in &RASSCF to the number of the deep core orbitals. Do you think so?

2. Which do you think are the deep core orbitals in the NUN(nitrogen-uranium-nitrogen) molecule, N([He]) and U([Xe]) or N([He]) and U([Kr])?

3. Without the “frozen” keyword in &RASSCF (in this case, frozen=0, all electrons minus twice the number of inactive orbitals is equal to the total number of active electrons), and without the “frozen” keyword in &CASPT2, and then run a single-point CASPT2 calculation. Will I find the default value of “frozen” in &CASPT2 in the ProjectName.log file?

Thanks a lot!

Best regards,
huangtf

Last edited by huangtf (2017-05-02 08:55:06)

Ignacio · 2017-05-02 14:33:55

huangtf wrote:

1. “The default is to freeze the max of those that were frozen in the RASSCF calculation and the deep core orbitals.” in the paragraph describing the meaning of the keyword “frozen” in &CASPT2 (page 128 in the manual) implies that, without the “frozen” keyword in &CASPT2, the default value of “frozen” in &CASPT2 will be equal to the set value of “frozen” in &RASSCF when the latter is equal to or more than the number of the deep core orbitals, and will be more than the latter when the latter is less than the number of the deep core orbitals. To keep that the default value of “frozen” in &CASPT2 is equal to the set value of “frozen” in &RASSCF and to avoid that freezing too many orbitals looks like I am freezing my whole system and thus &CASPT2 is actually doing nothing, I should set the value of “frozen” in &RASSCF to the number of the deep core orbitals. Do you think so?

No. You should have no frozen orbitals in &RASSCF and let &CASPT2 choose the default for you. Unless you have any reason to do otherwise, which doesn't seem to be the case.

2. Which do you think are the deep core orbitals in the NUN(nitrogen-uranium-nitrogen) molecule, N([He]) and U([Xe]) or N([He]) and U([Kr])?

Best way to find out is letting &CASPT2 select the frozen orbitals.

3. Without the “frozen” keyword in &RASSCF (in this case, frozen=0, all electrons minus twice the number of inactive orbitals is equal to the total number of active electrons), and without the “frozen” keyword in &CASPT2, and then run a single-point CASPT2 calculation. Will I find the default value of “frozen” in &CASPT2 in the ProjectName.log file?

Yes, for instance:

++    Orbital specifications:
      -----------------------
 
      Symmetry species                           1
                                                 a
      Frozen orbitals                            4
      Inactive orbitals                          6
      Active orbitals                           10
      Secondary orbitals                       128
      Deleted orbitals                           0
      Number of basis functions                148

huangtf · 2017-05-08 14:26:38

Dear Ignacio,

Thanks a lot for your kind and helpful reply! I have five questions to consult you.

1. I should have no frozen orbitals in &RASSCF, i.e. I don't write the “frozen” keyword in &RASSCF (by default, frozen=0). In this case, the number of inactive orbitals (the set value of “inactive” keyword in &RASSCF) equals all electrons minus the total number of active electrons, and divided by two. Do you think so?

2. According to the last sentence “With ANO type basis sets it is actually preferred to use the DK Hamiltonian and ANO-RCC in all your calculations.” in section 10.7.1 of the manual (page 511), it seems that scalar relativistic effects should be treated in my calculations by using the Douglas-Kroll-Hess (DKH) Hamiltonian (using the keyword “Douglas-Kroll” in &Seward). However, according to the sentences “The basis have been contracted using the Douglas-Kroll Hamiltonian and should therefore only be used in calculations where scalar relativistic effects are included. Seward will automatically recognize this and turn on the DK option when these basis sets are used.” interpreting the bold words “Relativistic ANO basis sets — ANO-RCC” in page 369, it is plausible that I don't need to include the “Douglas-Kroll” keyword in the &Seward input, because my basis set is ANO-RCC-VTZP. How should I do?

3. The keyword “RXXPyy” is employed to “Request arbitrary scalar relativistic Douglas-Kroll-Hess (DKH) correction to the one-electron Hamiltonian and the so-called picture-change correction to the property integrals (multipole moments and electronic potential related properties).” and “Recommended orders and parametrization is R02O02.” (see page 336). Do I need to use the “RXXPyy” keyword in &Seward?

4. By default, can local exchange (LK) screening be used automatically in &RASSCF and &CASPT2 in combination with Cholesky decomposition to reduce the cost of the two electron integrals? (only use the keyword “Cholesky” in &Seward in my input file)

5. Using the keyword “Cholesky” in &Seward means that, by default, using the Cholesky decomposition representation of electron repulsion integrals (two-electron integrals) in all stages of my calculations (integrals, CASSCF, and CASPT2), or that, applying the CD-CASPT2 approach in my calculations by default. Do you consider what I understand is correct?

Thank you in advance!

Sincerely yours,
huangtf

Ignacio · 2017-05-08 17:04:46

huangtf wrote:

1. I should have no frozen orbitals in &RASSCF, i.e. I don't write the “frozen” keyword in &RASSCF (by default, frozen=0). In this case, the number of inactive orbitals (the set value of “inactive” keyword in &RASSCF) equals all electrons minus the total number of active electrons, and divided by two. Do you think so?

Yes

it is plausible that I don't need to include the “Douglas-Kroll” keyword in the &Seward input, because my basis set is ANO-RCC-VTZP. How should I do?

Right, if you use ANO-RCC you don't need to add any special keyword to &SEWARD. Check the output if you want to make sure, it will say something like:

                  Relativistic Douglas-Kroll-Hess integrals:
                    - Parametrization         : EXP
                    - DKH order of Hamiltonian: 2
                    - DKH order of Properties : 0
                         - multipole moment operators
                         - electric potential operators
                         - contact operators
                  Atomic mean-field integrals

Do I need to use the “RXXPyy” keyword in &Seward?

No (see above)

4. By default, can local exchange (LK) screening be used automatically in &RASSCF and &CASPT2 in combination with Cholesky decomposition to reduce the cost of the two electron integrals? (only use the keyword “Cholesky” in &Seward in my input file)

There are many algorithmic tweaks for Cholesky. If you don't know/want to fiddle with them, just use "RICD" in &GATEWAY and forget the rest (i.e., no "Cholesky" in &SEWARD)

5. Using the keyword “Cholesky” in &Seward means that, by default, using the Cholesky decomposition representation of electron repulsion integrals (two-electron integrals) in all stages of my calculations (integrals, CASSCF, and CASPT2), or that, applying the CD-CASPT2 approach in my calculations by default. Do you consider what I understand is correct?

Yes, once you use CD (through "RICD" or "Cholesky") in &SEWARD, all later steps are forced to use CD, simply because conventional two-electron integrals are not available. You don't have to specify anything else.

huangtf · 2017-05-26 16:06:24

Dear Ignacio,

Thank you very much for your kind and helpful reply last time! I have three questions to consult you.

1. “Some care has to be taken in choosing the input orbitals, especially for the weakly occupied ones. Different choices may lead to convergence to different local minima.” (appearing in the sixth line of the second paragraph in page 294) Thus, I wonder, if the active space is “wrong”, $Project.GssOrb, $Project.ScfOrb, and $Project.RasOrb, which kind is the best starting orbitals used to choose the proper active space? Is it up to the total number of pairs of MOs that will be exchanged because of the specific chemical properties I want to investigate? (the less, the easier?)

2. “Always check that your active space did not change”, which means that, during the whole process of my CASPT2 geometry optimization, I should often inspect my active space, i.e. I can run the CASPT2 optimization and regularly open the $Project.rasscf.molden file while it is running (but after it's finished &RASSCF). When I find the active space become “wrong”, according to “In particular the type index can be changed directly in the file if the RASSCF program has converged to a solution with wrong orbitals in the active space.” (see p 63) I should kill the CASPT2 optimization job at once (generally, “kill PID number”), and then use the last $Project.RasOrb in which the ordering of MOs will be altered and the last geometry structure to do the CASPT2 optimization again. Do you consider that I comprehend it correctly?

3. I've written “>EXPORT MOLCAS_MOLDEN=ON” on my input file, and basis=ano-rcc-vtzp, but the following words appear in the $Project.log file:
“Sorry, Molden does not know how to handle
functions with angular momentum larger than g”
The $Project.GssOrb file is generated, the corresponding $Project.guessorb.molden file is not produced. (It's easy to speculate that the following $Project.rasscf.molden, $Project.geo.molden, and $Project.freq.molden will also not be outputted.) How to solve it, please?

Thank you for taking the time to read my problems! Please help me out.

With best wishes,
huangtf

Last edited by huangtf (2017-05-26 16:08:35)

Ignacio · 2017-05-27 10:47:04

huangtf wrote:

if the active space is “wrong”, $Project.GssOrb, $Project.ScfOrb, and $Project.RasOrb, which kind is the best starting orbitals used to choose the proper active space? Is it up to the total number of pairs of MOs that will be exchanged because of the specific chemical properties I want to investigate? (the less, the easier?)

Check these tips: http://www.molcas.org/documentation/man … 0000000000

Typically: pi and pi* orbitals, plus any sigma and sigma* orbitals if of bonds that break/form, and possibly lone pairs and d/f orbitals in metal centers.

When I find the active space become “wrong”, according to “In particular the type index can be changed directly in the file if the RASSCF program has converged to a solution with wrong orbitals in the active space.” (see p 63) I should kill the CASPT2 optimization job at once (generally, “kill PID number”), and then use the last $Project.RasOrb in which the ordering of MOs will be altered and the last geometry structure to do the CASPT2 optimization again. Do you consider that I comprehend it correctly?

Yes. You can either alter the orbitals or start from the scratch, but with the latest geometry in both cases, otherwise you gain nothing.

“Sorry, Molden does not know how to handle
functions with angular momentum larger than g”
The $Project.GssOrb file is generated, the corresponding $Project.guessorb.molden file is not produced. (It's easy to speculate that the following $Project.rasscf.molden, $Project.geo.molden, and $Project.freq.molden will also not be outputted.) How to solve it, please?[

The message is pretty clear. Molden does not support those basis functions, so you cannot use Molden with that system. You could use GRID_IT plus GV/LUSCUS. Or you could use a smaller basis that doesn't have h functions.

Molcas Forum

Announcement

#1 2017-04-25 09:10:46

Write input file for caspt2 optimization

#2 2017-04-25 09:56:45

Re: Write input file for caspt2 optimization

#3 2017-04-25 16:09:18

Re: Write input file for caspt2 optimization

#4 2017-04-26 08:24:19

Re: Write input file for caspt2 optimization

#5 2017-04-26 13:15:58

Re: Write input file for caspt2 optimization

#6 2017-04-26 17:20:15

Re: Write input file for caspt2 optimization

#7 2017-04-28 15:30:10

Re: Write input file for caspt2 optimization

#8 2017-04-28 16:00:34

Re: Write input file for caspt2 optimization

#9 2017-05-02 08:46:13

Re: Write input file for caspt2 optimization

#10 2017-05-02 14:33:55

Re: Write input file for caspt2 optimization

#11 2017-05-08 14:26:38

Re: Write input file for caspt2 optimization

#12 2017-05-08 17:04:46

Re: Write input file for caspt2 optimization

#13 2017-05-26 16:06:24

Re: Write input file for caspt2 optimization

#14 2017-05-27 10:47:04

Re: Write input file for caspt2 optimization

Board footer