Firefly and PC GAMESS-related discussion club

Learn how to ask questions correctly  
We are NATO-free zone

Re^8: MCSCF for a large system.


Dear Sanya,
Thanks for your detailed letter. Itís a big of misfortune that you are not a part of Dr Granovsky team Ė the ffly manual would be much better thenÖ
You stressed me much about wstate(1)=1,1 plan but havenít pursuit that wstate variation would not work. Of course I am not in position to argue against high-skilled Qchemist on this technical moment. Let it be my scientific religion at this point. What is the evident problem of your position is the lack of experimental data Ė are there any references about variation of wstate for the systems with constrained geometry? If there is nothing then you are the scientific follower either, yet the version of your religion is different :). The good news is that scientific religions are verifiable though. So I have a proposition: why not to setup two approaches and see what modeling technique gives a better result in experiments with real compounds and protein. In the end youíll have experimental fact about wstate variation with one of the most adequate exotic system and Iíll got some compounds if lucky enough. The only things I need from you is qualified advises about wstate variation approach (the win would give you nothing if you know that I did many stupid mistakes through the modeling). So please imagine the best computational experiment for wstate variation and briefly describe it; my obligation will then be to report about the results (computational and experimental). Let the experiment judge us! Does it sound like a plan?
If you are interested in it, the current protocol is:
Branch A, internal:
RHF B3LYP optimization of a big protein fragment then cutting and freezing of atoms of the heart of active site
following GVB and following MCSCF geometry optimization with constrained all atoms but those that participate reaction (reaction coordinate is constrained and sampled with 0.1A step)
for mcscf part, the geometry optimization is made on a grid (wstate vs coordinate of reaction), grid is to be of manually tractable size: 2 states, wstate ratio is 0.1, 0.3, 0.5, 0.7, 0.9 (to fit golden ratio and halves on a regular grid).
Branch B, entropic (rather FYI, I am a sort of expert here):
mcscf parameterization of reagents and (intermediate) products complexes from Branch_A for MM MD application (GROMOS/GROMACS here) (normal modes, torsions fit, electric moments fit)
measure protein free energy with alchemical slow-growing perturbation method in MD and extract reaction coordinate distribution over perturbation parameter.
Joining of Branch_A and Branch_B (probably require more thoughts):
update dH/dlambda in Branch_B (protein free energy) with data from QM Branch_A with finite differences of mcscf energies at different reaction coordinates
consider how Branch_A and Branch_B fits the experiment alone and (several) modes of their mixtures.
P.S. You of course realize that it takes some time (several months) to carry out computations and we are limited to c.a. 350K existing compounds from our partners (itís c.a. 30 items only to screen&test considering geometry of the protein cavity and electronic properties of the substances).

On Tue Sep 30 '14 4:04pm, sanya wrote
>First of all, I should say: orbitals are not intrinsic properties of molecules. This is why they cannot be taken as molecular descriptors for sceening. Change the basis set, change the functional -- and you get different orbitals. The electron density is the intrinsic property. What is more, it is observable. However, you cannot predict a priori how the observed density of the molecule changes upon electron loss or capture. Therefore, you cannot predict the redox behavior of the molecule without simulating the target reaction in more or less realistic conditions.

>Statements like "The HOMO of molecule A has more contribution from the N in 3-position than the HOMO in molecule B; therefore, A is more active in the reaction" can make some sense only when (1) single-reference method like HF or DFT is used, (2) the basis set and functional (for DFT) are the same for A and B, (3) the nodal structure of HOMO in A and B is the same, which may not be the case even if A and B are very similar.

>>Should one expect significant orbital perturbation due to transfer of beta electron only in a system of molecules with constrained geometry? I donít think so.

>Electron transfer itself _is_ a substantial orbital perturbation. And geometry constraints are much less than you think: below, you say that O-O bond breaks, and this is not the only change, I'm sure. The electron density change (from the density of N electrons on M orbitals you get N+1 electrons on M+1 orbitals, which is definitely not the same :) ) changes the potential surrounding the nuclei and causes them to adopt new geometry. Although the geometry changes may not include large conformational changes, you'll see them, for example, as the changes in the interatomic distances and bond lengths, even in constrained environments.

>>What would change thus (as the reaction proceeds)? The configurations rate I believe. Does it make sense now?

>Of course, not.
>You must be confusing orbital configurations and states. Let's make it clear.
>First of all, we have an atomic basis set (atomic orbitals, AOs), that is, the set of functions (more or less resembling hydrogen-like orbitals of free atoms) centered on each atom of the molecule.
>In HF or DFT the wavefunction (or electron density) is described by the linear combinations of these atomic basis functions. These linear combinations are called molecular orbitals (MOs), and the wavefunction (or its formal analog in DFT) is constructed as Slater determinant. The coefficients of these linear combinations (contributions of each AO) are chosen so as to minimize the matrix element of the Hamiltonian with the corresponding Slater determinant. The total electron density, which is the observable, is governed by the _occupied_ orbitals. Therefore, adding or withdrawing an electron _does_ change the electron density.

>That is, at the first step, we have a Slater determinant built of molecular orbitals, each of them being a linear combination of atomic orbitals (sort of The House That Jack Built ;) ). In many cases, the ground state of a molecule can be well described by a single Slater determinant.

>If we take all (not only occupied) molecular orbitals obtained at the previous step (usually, HF orbitals are taken) and construct all possible Slater determinants in which one or more occupied orbitals are replaced by the same amount of unoccupied ones, this is interpreted as excitations (we excite electrons from occupied orbitals to unoccupied ones). However, the excited states cannot be described by individual Slater determinants (configurations). Each excited state is described by a linear combination of excited Slater determinants, and the coefficients of this combinations are again chosen so as to minimize the matrix element of the Hamiltonian. This is one more floor in the House That Jack Built. An excited state can have one dominant configuration, but not necessarily. This method is called Configuration Interaction (CI). It can be of different order (CIS, CISD, CISDT, etc.) or include all possible excitations within a given set of occupied and unoccupied orbitals (active space). The latter version is called Complete Active Space CI, or CAS-CI. If the active space includes all occupied and unoccupied orbitals, corresponding CAS-CI becomes Full CI.

>Multiconfigurational SCF includes adjustment of molecular orbitals in addition to the adjustment of CI coefficients. The larger the active space (i.e., the closer the active space to the full space), the less orbital adjustment is needed. And orbital adjustment is nothing but changes in AO coefficients in MOs. This is the attic of old Jack's House ;) Therefore, the only unchanged thing is basis set.

>In principle, we could adjust the orbitals for each state individually, and these orbital sets would be different for each state. But the only criterion for adjustment of MO and CI coefficients is minimization of the matrix element of the Hamiltonian, which is valid only for the ground state of the given multiplicity and symmetry. Therefore, if we need an excited state of the same mult. and symmetry as the ground state, we cannot adjust its wavefunction individually. We have to use state-averaging.

>State-averaging is the way to obtain the orbitals equally good (or equally bad) for several states. In this case, the Hamiltonian is built with the state-averaged density matrix rather than with individual MOs, or configurations, or states. And the value to minimize is the average of the Hamiltonian with this density matrix. The weights of the states here are _not_ adjusted (not variate) and, therefore, should not be interpreted. But if you change the states being averaged, you'll definitely get different state energies and state wavefunctions. And the wavefunction in this case is the linear combination of Slater determinants build on MOs, which are, in their turn, linear combinations of AOs.

>>The experimental part of my project is screening of cofactor candidates to speedup the reaction. Thus orbital things go into play much later, after cofactor initiated the oxygen (surely I see the shifts as breaking O-O bond, but they are not the knowledge I need to predict cofactor efficiency).

>Careful study of the reaction mechanism with MCSCF will help you find adequate descriptors for screening.

[ Previous ] [ Next ] [ Index ]           Wed Oct 1 '14 8:47pm
[ Reply ] [ Edit ] [ Delete ]           This message read 660 times