Ceilidh ... Re^9: MCSCF for a large system.

Firefly and PC GAMESS-related discussion club

Learn how to ask questions correctly

We are NATO-free zone

Re^9: MCSCF for a large system.

Dear Sanya,
Your approach to manual writing is very extraordinary: give users to get lost in the ugly documentation first and then fix the places where ones stack the most. It’s the approach if one wona to melt stuff into pure diehard qchemists and a slag but would be a bit repulsive to an average novices.
Understanding that you had passed through the above ‘metallurgy’, I am emphasizing the need for a good scientist to stay open-minded (either it is religion in science or wstates balance problem). It’s of course up to you: join the following discussion or mark the problem as a wasting but I am screening the transition of reaction with two wstate variation followed by geometry optimization. The results are expected by Mon/Tue, it will a separate topic. Why not to check how the variation looks like in energy terms of a real problem, are the routes through wstates make any difference or it’s just a matter of mkcals etc. In perspective, with geometry statistics from “a MD on the PESs parameterized by QM” we will see what states combination the Lord (or just the protein :) ) prefers at each reaction coordinate.
P.S. Thanks for the suggestion about reaction force field. Its however maybe not the best way to proceed. We don’t need a field for model a chemical reactions as already making the FF which is tuned to accurately model only this particular reaction (this what the overwhelming citation above stands for, isn’t it?). Of what we really need is a good exploration of protein configuration space. Thus GROMOS which is about protein-solvent equilibriums would be superior over FFs focused on reaction description.

Regards,
Alex

On Thu Oct 2 '14 0:54am, sanya wrote
------------------------------------
>>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…

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

Fri Oct 3 '14 9:41pm

This message read 587 times