Firefly and PC GAMESS-related discussion club

 
Learn how to ask questions correctly  
 
 

 

sanya
sanya@photonics.ru

Hmm, actually, I am a member of Dr Granovsky team :) and it's my fault that I didn't pay enough attention to MCSCF part of the manual. I think, we'll include my explanation into the manual after some editing.

My experience shows that program manuals are the best quantum chemistry textbooks. However, they may include some practical hints, but completely lack rigorous theory. Therefore, manuals cannot be a substitute to fundamental textbooks like those of, say, A. Szabo & N. Ostlund, T. Helgaker et al., O. Sinanoglu, etc. There are some good books in Russian by V.I. Pupyshev & N.F. Stepanov, V.G. Tsirel'son, etc.

>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.

Religion is something opposed to science ;) Every statement in science should be either rigorously proven or be based on experience rather than taken on trust.

I state that wstate is not a variational parameter (unlike MO and CI coefficients), because, once chosen, it cannot be optimized with respect to any reasonable criterion. It is not an experimental parameter either, because there is no physical property that can be directly related to wstate. It is just a helpful way to obtain orbitals describing several states on equal grounds. If you don't like wstate(1)=1,1, you may take wstate(1)=1,0 or 1,0.5, or whatever you like. But you cannot change wstate along the reaction path, because it would be equivalent to, say, changing DFT functional from BLYP to B3LYP to BHHLYP along the reaction path. The results will be meaningless.
All the above (1) follows from the mathematical formulation of SA-CASSCF method, (2) was tested by experiments many times.

>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?

Of course, not. Because, to get the idea of what you're going to do, you may first try to calculate the reagents by BLYP, the transition state by B3LYP, and the products by BHHLYP, and vice versa, and try some other combinations of three functionals with three geometries. Then, compare the resulting reaction profiles. This will demonstrate you all the inadequacy of this approach much faster.

>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)

This seems OK

>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).

From what I read above, I see that you still haven't digested my explanation about orbitals, configurations, and states. Probably, it was not clear enough :(

You don't have to change wstate. Instead, you can plot, for example, the weights of certain configurations in the states of interest vs. reaction coordinate. If you are assiduous enough, you can even plot the weights of certain AOs in the MOs whose population changes most noticeably during the reaction. The weights of configurations and the weights of AOs are variational parameters, and they can be related (indirectly, unfortunately) to the observed electron density.
On the contrary, wstate is a sort of a voluntary parameter. If you need two states equally, you should give wstate(1)=1,1. Giving preference to one of the states means only that you like this state more than the other one. That is, if you set wstate(1)=0.1,0.9 for the reagent geometry and wstate(1)=0.9,0.1 for the product geometry, this only means that you like excited state of the reagent more than its ground state and you like the ground state of the product more than its excited state.
But one more important difference between the science and religion is that in science _your_attitude_does_not_matter_. You cannot like one state and hate the other. You can only study them.

And why do you talking all the time about constrained geometry of the active site? If you properly perform _relaxed_ scan along the reaction path, you'll be surprised by the large reorganization. Don't forget, relaxed scan rather than simple scan!

>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.

So, you're going to perform a MD on the PESs parameterized by QM. This seems sound. Moreover, you're going to simulate surface hopping from one PES to the other, because you're expecting that you PESs cross in some region. This is difficult but possible. I never tried this, but there are some good papers (sorry, I cannot send them to you right now, but you may search Web of Science and get them faster than I find them in my computer).

The problem with the parameterization is that the reaction coordinate involves not the torsions, but the bond lengths. Therefore, you'll need a reactive force field. I'm sure you're fully aware of this.

------

Ufff, this explanation made me feel as if I've turned myself inside out...