Dear friends,
Recently, I begin to follow the geometry optimization tutorials (https://github.com/psi4/psi4numpy/blob/master/Tutorials/13_Geometry_Optimization/13a_Internal-Coordinates-Bmatrix.ipynb) of the psi4numpy. In this tutorial, the Ammonium molecule has been chosen as the example. In the tutorial, the author gave the codes about how to calculate redundant coordinate, such as bond, angle. However, for general, dihedral torsion should be calculated. Unfortunately, there is no codes about how to calculate it. Can you provide the general codes about calculating redundant coordinates from cartesian coordinates for me?
Yes, I have read carefully this function in opt_helper. I can understanding it and know how to calculate it. But I don’t know how to judge four atoms to form a dihedral angle.
I don’t know what it means to “judge four atoms to form a dihedral angle”. Could you try saying it a different way? Do you mean “how to know when the dihedral angle formed by four atoms is a good internal coordinate”?
Dihedral angles that are fixed by symmetry. For formaldehyde, D(H1,C,O,H2) must be 180 degrees.
Dihedral angles for consecutively bonded atoms. For sulfuric acid, H-O-S=O would be a good choice because H is bonded to O, O is bonded to S, and S is (double) bonded to that other O.
Dear jmisiewicz,
According to your suggestions, I have writen some codes to calculate the “good dihedrals”. Unfortunately, I met some trouble (I have created successfully bond, angle). For example, I want to create redundant coordinate of the “methane”. According to your suggestion 2, there is no consecutively four bonded atoms. The bond matrix is:
However, in fact, there exists 4 dihedrals (extracted from Gaussian program output), i.e., D(2,1,4,3), D(2,1,5,3), D(2,1,5,4), D(3,1,5,4).
I hope that you can use a more general molecule as the case in your tutorial of psi4numpy.
Psi4 is so good as Python interface. Psi4Numpy is a good reference for studying the basic quantum theory. I learn much knowledge following your tutorial of Psi4Numpy.
Automatically determining appropriate internal coordinates for a general molecule is an extremely difficult problem. The considerations I gave are appropriate for most molecules, but not all. If you really want more information on this, ask @Rollin_King.
Taiping, you can definitely add torsions if you wish manually using analogous function calls, for example,
intcos.append( tors.TORS(i, j, k, l) )
It is just unclear, as the other folks responded, why you would want to do this. If there is no bonding pattern i-j-k-l, then the coordinate is not effective in geometry optimizations (because the Hessian/2nd derivative is very not diagonal in such coordinates). This highlights one of the problems with Z-matrices, which force the use of dihedrals.
Dear Rollin_King,
Thank you for your response.
In fact, I want to calculate Wilson B matrix according to the cartesian coordinates(the current geometry is optimized by other quantum program). Firstly, therefore, I want to calculate redundant coordinates and the transformation matrix between redundant and cartesian coordinates, i.e., Wilson B matrix.
These are things you may already know, but I’ll point out in case you do not know these:
In order to define a Wilson B matrix, you need a geometry in Cartesians and also a choice of internal coordinates.
The Wilson B matrix is technically not a transformation matrix between redundant and cartesian coordinates. It is the best linear approximation to the nonlinear map between Cartesian displacements and internal displacements. It is a map between changes in coordinates, not coordinates themselves, and it is the matrix approximation to a much more complicated function.
Also, why do you need a Wilson B matrix in the first place? The large problem you are trying to solve may have a simpler solution.
Dear jmisiewicz,
I want to calculate the Duschinsky matrix between two sets of normal modes. I will read force constant matrix and frequency information from other quantum chemical software. The final purpose are to calculate Duschinsky matrix and decomposite reorganization energy. The reference is “J. Chem. Phys. 115, 9103 (2001)”. Can you provide some suggestions for me? Thanks a lot.
I’m not familiar with these matrices, but I’ve skimmed the paper.
You seem to be trying to create a Python program to take in Cartesian force constant data for an arbitrary molecule at two different stationary points and compute the B matrix, as an intermediate for computing Duschinsky matrix. Am I understanding you? We are quickly getting to the limit of what I know, but I will have different things to say if you want to compute for a particular molecule or to create a program to do this for general molecules.
Automatically generating appropriate non-redundant internal coordinates is a hard problem, and one I’m not qualified to help with. I see three options:
Have the user manually create the internal coordinates.
In cases like methane where there are not enough “good” internal coordinates, create ones that aren’t “good.” Create as many as you need, and then use the procedure from Sec. IIA of the paper to remove redundancies.
Abandon internal coordinates and construct the matrix in mass-weighted Cartesians. This may or may not be possible. Again, I’m not familiar with Duschinsky matrices.