# Gaussian 16: Relaxed scan using Jacobi coordinates expressed using generalized internal coordinates

For teaching porpouses, I would like to make a relaxed scan for an A-BC system (for example, the isomerization of of H-NC to H-CN or CH3-NC to CH3-CN) with Gaussian 16 using generalized internal coordinates. In particular, the coordinates to be scanned would be the distance of atom A to the center of mass of the BC, R, and θ, the angle between R and the direction of the BC bond (rAB) as illustrated in the figure below.

Could anybody give me an example of how to define these particular coordinates?

Following the suggestions of @S.R.Maiti I have prepared an input for H-NC to scan the R coordinate at a fixed θ angle, and scanning the θ angle or both coordinates. The results demonstrate a mixed successes.

Sample code for the can of the R coordinate

This kind of scan completes withour error.



Test GIC HCN

0 1
C -0.6174 0.0000 0.0000
N  0.5290 0.0000 0.0000
H  1.1954 0.1046 0.0000

comX = (X(1)*12.0 + X(2)*14.0)/(12.0+14.0)
comY = (Y(1)*12.0 + Y(2)*14.0)/(12.0+14.0)
comZ = (Z(1)*12.0 + Z(2)*14.0)/(12.0+14.0)
comToBX = X(2) - comX
comToBY = Y(2) - comY
comToBZ = Z(2) - comZ
ComToBdist = SQRT(comToBX**2 + comToBY**2 + comToBZ**2)
comToCX = X(3) - comX
comToCY = Y(3) - comY
comToCZ = Z(3) - comZ
ComToCdist = SQRT(comToCX**2 + comToCY**2 + comToCZ**2)
dot = comToBX*comToCX + comToBY*comToCY + comToBZ*comToCZ
BcomC(Freeze) = ARCCOS(dot / (ComToBdist * ComToCdist))
ComToCdist(NSteps=13,StepSize=0.092)


Sample code for the scan of the R and θ coordinates

When I tried to scan the θ coordinate (on its own or in a double scan), I got an error message, I guees, related to the definition of the angle (see below).



Test GIC HCN

0 1
C -0.6174 0.0000 0.0000
N  0.5290 0.0000 0.0000
H  1.1954 0.1046 0.0000

comX = (X(1)*12.0 + X(2)*14.0)/(12.0+14.0)
comY = (Y(1)*12.0 + Y(2)*14.0)/(12.0+14.0)
comZ = (Z(1)*12.0 + Z(2)*14.0)/(12.0+14.0)
comToBX = X(2) - comX
comToBY = Y(2) - comY
comToBZ = Z(2) - comZ
ComToBdist = SQRT(comToBX**2 + comToBY**2 + comToBZ**2)
comToCX = X(3) - comX
comToCY = Y(3) - comY
comToCZ = Z(3) - comZ
ComToCdist = SQRT(comToCX**2 + comToCY**2 + comToCZ**2)
dot = comToBX*comToCX + comToBY*comToCY + comToBZ*comToCZ
BcomC = ARCCOS(dot / (ComToBdist * ComToCdist))
ComToCdist(NSteps=13,StepSize=0.092)
BcomC(NSteps=13,StepSize=5.0)


Once Gaussian has completed the scan for the first value of θ, it produces the following error message:

                          !   Optimized Parameters   !
! (Angstroms and Degrees)  !
--------------------------                            --------------------------
! Name       Definition         Value          Derivative Info.                !
--------------------------------------------------------------------------------
! comX       GIC-1             -0.3944         -DE/DX =    0.019               !
! comY       GIC-2             -0.0342         -DE/DX =    0.0038              !
! comZ       GIC-3              0.0            -DE/DX =    0.0                 !
! comToBX    GIC-4              1.0547         -DE/DX =    0.0124              !
! comToBY    GIC-5             -0.0005         -DE/DX =    0.0017              !
! comToBZ    GIC-6              0.0            -DE/DX =    0.0                 !
! ComToBdist GIC-7              1.0547         -DE/DX =    0.0124              !
! comToCX    GIC-8              3.4508         -DE/DX =   -0.0482              !
! comToCY    GIC-9              0.3002         -DE/DX =   -0.0096              !
! comToCZ    GIC-10             0.0            -DE/DX =    0.0                 !
! ComToCdist GIC-11             3.4638         -DE/DX =   -0.0488              !
! dot        GIC-12             3.6392         -DE/DX =   -0.0074              !
! BcomC      GIC-13             0.0873         -DE/DX =   -0.0032              !
--------------------------------------------------------------------------------
NOTE: GIC-type coordinates are in arbitrary units.
Iteration  1 RMS(Cart)=  0.03691776 RMS(Int)=  1.36059875...
Iteration100 RMS(Cart)=  0.00247220 RMS(Int)=  1.49252761
New curvilinear step not converged.
Error imposing constraints
Error termination via Lnk1e in C:\G16W\l103.exe at Sat Dec 09 22:17:54 2023.
Job cpu time:       0 days  0 hours  4 minutes 39.0 seconds.
Elapsed time:       0 days  0 hours  2 minutes 34.5 seconds.
File lengths (MBytes):  RWF=     11 Int=      0 D2E=      0 Chk=      6 Scr=    1


I guess the problem is related to the definition of the angle. Is there a way to avoid it?

• Hmm, only the ComToCdist and BcomC should be GIC's I think? The rest should be the usual primitives (bonds, angle etc.). I haven't used Gaussian in a while, so I don't know what is the right keyword to do this. The other variables (comX, comY etc.) are just needed to define the needed angle and COM, they should not be GICs Commented Dec 9, 2023 at 22:16
• Also, the BcomC angle is 0.08 radian from the printout, which means that the B matrix element of this angle will go to infinity, so the internal to cartesian backtransform will be difficult to converge. (This is a problem with the arccos definition, the derivative of angle has sin theta at the denominator, and sin theta -> 0 as the angle goes to zero) Commented Dec 9, 2023 at 22:21

Based on the GIC page in Gaussian website and your picture, the definition of the coordinates would be:

comX = (X(A)*mA + X(B)*mB)/(mA+mB)  # these are the cartesian components
comY = (Y(A)*mA + Y(B)*mB)/(mA+mB)  # of the centre of mass of A-B system
comZ = (Z(A)*mA + Z(B)*mB)/(mA+mB)

# for the angle
comToBX = X(B) - comX
comToBY = Y(B) - comY
comToBZ = Z(B) - comZ
ComToBdist = SQRT(comToBX**2 + comToBY**2 + comToBZ**2)
comToCX = X(C) - comX
comToCY = Y(C) - comY
comToCZ = Z(C) - comZ
# below is the COM --  C distance coordinate as defined in the picture
ComToCdist = SQRT(comToCX**2 + comToCY**2 + comToCZ**2)

dot = comToBX*comToCX + comToBY*comToCY + comToBZ*comToCZ
# below is the B - COM - C angle defined in your picture
BcomC = ARCCOS(dot / (ComToBdist * ComToCdist))


You have to replace the following variables in the above expressions with numbers (constants): mA = mass of atom A, mB = mass of atom B, A = index of atom A, B = index of atom B, C = index of atom C. (also remove the comments)

Note that the scanning using these coordinates will fail when the angle $$\theta$$ (B-COM-C) is close to 0 or 180 degrees. This is due to the B matrix elements of the corresponding angle going to infinity as the angle approaches linear or zero angle. This is a well-known problem with the angles defined with arccosine of the normalised dot product (which is the way I have written it here). It's quite difficult to solve this frankly... most softwares replace the angle coordinate with a linear bend coordinate (that requires a dummy atom or a dummy axis) when the angle is close to linear, but that would be counterproductive here since you want to scan the whole range of angle.

• Thanks @S.R.Maiti. Working on it. I had the vain hope that there was a keyword in Gaussian that provided the com coordinates without calculations.
– PAEP
Commented Dec 3, 2023 at 18:28

Defining the geometry in terms of z-matrix internal coordinates would make this very easy, I think. Define COM as a dummy atom, place another dummy atom to define the molecular plane, place A and B with suitable distances to dummy-1 and 90 degree angles to dummy-2, and B having a torsional angle of 180 relative to A. Now define C with a distance to dummy-1 and an angle to B and torsional of 180 to A. Simply keep the first two distances for A and B fixed, and give the scan range and step size for the two variables defining C, as given in the Gaussian manual for the keyword SCAN. That should work for all combinations of R and theta, including 0 and 180.