13
$\begingroup$

I would like to perform a relaxed scan in Gaussian09 on my molecule using the modredundant optimization. Because I want to scan over two bond lengths simultaneously, I defined some dummy atoms which are close to the real atoms whose bond lenghts should be scanned. However, this approach does not work, because I can't properly index the dummy atoms. The Gaussian09 documentation says (emphasis mine):

[Type] N1 [N2 [N3 [N4]]] S nsteps stepsize [[min] max]]

N1, N2, N3 and N4 are atom numbers or wildcards (discussed below). Atom numbering begins at 1, and any dummy atoms are not counted.

How can I then perform a scan referencing a dummy atom? I don't want to do the scan using the z-matrix optimization because it does not converge.

Here is my input file:

%nprocshared=16
%mem=5gb
%chk=4MeAminoButanal
#P HF/6-31G(d) opt=(modredundant) scrf=(pcm) nosymm

4MeAminoButanal dummy scan

0 1
 C     0.000000     0.000000     0.000000
 N     0.000000     0.000000     1.451901
 C     1.325038     0.000000     2.055476
 C     1.991320     1.377897     2.032970
 C     1.229182     2.449098     2.809397
 C     1.309917     2.335179     4.306919
 X    -0.248373    -0.404008     1.611583
 X    -1.041414    -0.131882     1.066579
 H     1.218960    -0.323776     3.083767
 H     1.990264    -0.714465     1.564774
 H     2.993872     1.281489     2.436811
 H     2.106328     1.711600     1.006802
 H     1.602469     3.443240     2.566292
 H     0.178120     2.455789     2.536445
 H     0.708887     3.065770     4.854632
 H    -1.016528    -0.113770    -0.358506
 H     0.375985     0.943036    -0.379032
 H     0.606031    -0.798850    -0.431982
 X     1.641992     1.939406     4.607116
 X     1.311963     2.318294     5.471711
 O     1.974063     1.543637     4.907309
 H    -0.496745    -0.808017     1.771266

B 7 19 S 60 -0.05
B 22 7 F
B 2 7 F
B 6 19 F
B 21 19 F

When I run it, Gaussian09 obviously complains with

The following ModRedundant input section has been read: Invalid atom number. B 7 19 S 60 -0.05

Improve this question
$\endgroup$
6
  • 5
    $\begingroup$ This question might be better suited at Computational Science, but it is not off-topic here. You can indeed not use dummy atoms with a relaxed scan via modredundant, because they are not existent in that coordinate system. There is also no option, where this would be possible. Can you explain, what you are actually trying to do. To me it appears, that you are moving the middles of the CO and NH bonds closer together. If your input would work, then it would very likely not produce what you intended. (Btw. You are using too much proc and mem for HF.) $\endgroup$
    – Martin - マーチン ♦
    Commented Apr 14, 2016 at 11:36
  • 3
    $\begingroup$ @Martin-マーチン I do not agree with you: this problem is very specific for a QChem program, so it is very unlikely that Computational Science people can say much about it, while many Chemist have experience with the given software. However I fully agree with you that OP may want to do something where dummy atoms are most probably not needed at all. $\endgroup$
    – Greg
    Commented Apr 14, 2016 at 13:13
  • $\begingroup$ @Martin-マーチン That is exactly what I'm trying to do. How would you do it? $\endgroup$
    – shrx
    Commented Apr 14, 2016 at 14:04
  • $\begingroup$ @Greg From our tour: Don't ask about... Computational questions ... . Other than a technical aspect of the program Gaussian which is used to calculate molecules it does not have any bond to chemistry. Computational Science has the appropriate tag and most of the program users that are active here are (somewhat) active there. If you want to talk about this, we should probably go to meta. $\endgroup$
    – Martin - マーチン ♦
    Commented Apr 14, 2016 at 14:07
  • $\begingroup$ @shrx I actually meant what you want to investigate. What is the chemistry behind this approach. How I would go about it, you can see in my answer. I guessed what you want to do though. $\endgroup$
    – Martin - マーチン ♦
    Commented Apr 14, 2016 at 14:11

2 Answers 2

Reset to default
16
$\begingroup$

Disclaimer: This answer concerns the use of Gaussian 09. While the ModRedundant interface is retained and therefore the same restrictions also apply in newer versions, there is now a newer way of using internal coordinates. See my other answer for this option.

Like I posted earlier in the comments, you can indeed not use dummy atoms with a relaxed scan via the ModRedundant interface, because they are not existent in that coordinate system. Internal redundant coordinates are generated only from actual atoms. I tried to find an explanation of that, but failed. I believe the reason for that behaviour is, that dummy atoms do not have basis functions and are therefore not included in the wave function.

In this case you need alternative means to what you actually want to achieve. That of course depends on what you would like to model. From the input you provided, I can see, that you are trying to move the middle points of the carbon-oxygen and nitrogen-hydrogen bonds closer together, possibly to simulate the carbon nitrogen bond formation and proton transfer. While I seriously doubt this will happen in one step, this is at least one way to try to go about it.

If this method had worked, I believe you would have gotten a different than expected result. There are some more issues, that you might have not considered with your approach.
You might think that the following two lines keep the nitrogen hydrogen bond distance fixed. They do not.

B 22 7 F
B 2 7 F

Because the bond angle A 22 7 2 is still flexible, and the bond length B 22 2 is still flexible, this can possibly do whatever it wants. The same does hold for the carbon oxygen bond.

B 6 19 F
B 21 19 F

If I am correct with my assumption, then what you want look at is the carbon nitrogen bond formation. Why not take this internal coordinate and scan it? It is easy and it should pretty much do what you like. The following input should take you there. Note that I changed the method to pm6 to massively cut down computer time - and even with that, it took to long for my patience. (I did also not include any trailing lines about resources - see my notes below.)

#P pm6
opt=(modredundant)
scrf=(pcm)
nosymm

4MeAminoButanal scan

0 1
C        0.000000     0.000000     0.000000
N        0.000000     0.000000     1.451901
C        1.325038     0.000000     2.055476
C        1.991320     1.377897     2.032970
C        1.229182     2.449098     2.809397
C        1.309917     2.335179     4.306919
H        1.218960    -0.323776     3.083767
H        1.990264    -0.714465     1.564774
H        2.993872     1.281489     2.436811
H        2.106328     1.711600     1.006802
H        1.602469     3.443240     2.566292
H        0.178120     2.455789     2.536445
H        0.708887     3.065770     4.854632
H       -1.016528    -0.113770    -0.358506
H        0.375985     0.943036    -0.379032
H        0.606031    -0.798850    -0.431982
O        1.974063     1.543637     4.907309
H       -0.496745    -0.808017     1.771266

B 6 17 F
B 2 18 F
B 2  6 S 60 -0.05

Also note, that I did not check that the N-C bond length will become shorter than an Ångström, so the last steps take forever to converge and it will die in the last step running out of cycles. Run it with something more reasonable. The calculation itself produced a somewhat reasonable result. The following animation cycles through all optimised structures. It was made with ChemCraft.

scan animation

Some more notes on your input:

  • %nprocshared=16 is quite a number you are requesting for a Hartree-Fock calculation. Consider that parallelisation is not linear. While it became better throughout the years, I guess 4 threads is the most reasonable for this system size. You can go higher, I'd not recommend to go above ten though.
  • %mem=5gb is probably too much. I doubt that any of these processes will exceed 200MB of memory, probably a lot less.
  • %chk=4MeAminoButanal You should consider adding the .chk ending in all cases. Gaussian is known to mess things up when there are non-ordinary characters like _-. in a file name.
  • It is maybe a good idea to not write the whole route section in one line. Consider this as improving the readability of your input. You can cluster similar commands (like pop statements) in one line, etc..
    In previous versions, Gaussian might have had trouble reading in anything beyond the 80th char. I cannot reproduce this bug anymore and it's been a long time since I was told about it, so this should not be an issue. However, a well structured input file can help you finding possible syntax errors that result from typos.
Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ Gaussian sometimes has (still) trouble reading in anything beyond the 80th char. Wow, I never knew this... Can you offer up an explicit example? I'd love to see this behavior. Great answer by the way. $\endgroup$
    – LordStryker
    Commented Apr 14, 2016 at 14:42
  • 2
    $\begingroup$ @LordStryker Thank you and happy to see you. I could not reproduce this bug anymore, so I am fairly certain they fixed it. Gaussian used a fixed width of 80 char in the input (and still uses that as output for readability, but not everywhere) and break longer lines, apparently that could have destroyed some keywords. My former supervisor told and showed us, but that might have well been the 98 or 94 version. I'll fix that statement. $\endgroup$
    – Martin - マーチン ♦
    Commented Apr 15, 2016 at 5:07
  • $\begingroup$ I just realised, that this bug may have resulted from preprocessing of the input file, inserting spaces within a keyword, which the produces a syntax error. It might not be gaussians fault after all... $\endgroup$
    – Martin - マーチン ♦
    Commented Apr 15, 2016 at 5:26
  • 1
    $\begingroup$ Indeed I'm trying to do the C-N bond formation and simultaneous proton transfer. I tried doing it the simple way like you did, to just move the C and N atoms closer together. However, this does not lead to the desired proton transfer, as evident in your animation. I provided my input file as an example, but I am interested in a general way to do simultaneous scan over more than one variable with the modredundant optimization in Gaussian. For now, I have resorted to a not yet fully-functional python script that simulates a multi-variable scan. I can provide the code if there is interest. $\endgroup$
    – shrx
    Commented Apr 15, 2016 at 7:51
  • $\begingroup$ @shrx Well technically you are scanning only one parameter in your example input. You can do a multidimensional scan by providing more than one scan command in the modredundant section. I am still not sure what you are actually trying to achieve, because I believe that there is no true concerted transition state. $\endgroup$
    – Martin - マーチン ♦
    Commented Apr 15, 2016 at 8:14
8
$\begingroup$

Disclaimer: This answer concerns the use of Gaussian 16. While the ModRedundant interface is retained and therefore the same restrictions also apply to it as mentioned in the other answer. However, there is now a newer way of using internal coordinates as explained below.

The new magic is provided by the generalised internal coordinates (GIC), which provide awesome new ways of defining input. It is invoked via the geom or opt interface and the option specified is GIC. That alone doesn't let you perform scans, it just changes the internal use of coordinates and builds them. It has to be combined with the usual ModRedundant, or the new AddGIC (GIC would then be redundant) options.

The latter invokes the reading of a new input section after the molecule definition, just like in Gaussian 09. However, if you read through the info page, you will quickly notice, that many more options are possible, including mathematical operations. In this way you can define intermediate anchors which you may use as substitutes for dummy atoms. Dummy atoms in the geometry definition still don't count, and are completely superfluous with this approach.

The following input will perform the scan that you initially intended. It will break close to the end, as a four-membered ring will form and with two sides and the mid-points distance fixed, it simply runs out of cycles. You may want to adopt the input according to your needs. I'll explain a little more below the input-file.

#P PM6
OPT(MaxCycle=100,Loose)
SYMMETRY(none)
GEOM(ModRedundant,GIC)

4MeAminoButanal scan

0 1
C        0.000000     0.000000     0.000000
N        0.000000     0.000000     1.451901
C        1.325038     0.000000     2.055476
C        1.991320     1.377897     2.032970
C        1.229182     2.449098     2.809397
C        1.309917     2.335179     4.306919
H        1.218960    -0.323776     3.083767
H        1.990264    -0.714465     1.564774
H        2.993872     1.281489     2.436811
H        2.106328     1.711600     1.006802
H        1.602469     3.443240     2.566292
H        0.178120     2.455789     2.536445
H        0.708887     3.065770     4.854632
H       -1.016528    -0.113770    -0.358506
H        0.375985     0.943036    -0.379032
H        0.606031    -0.798850    -0.431982
O        1.974063     1.543637     4.907309
H       -0.496745    -0.808017     1.771266

xCO(inactive)=XCntr(6,17)
yCO(inactive)=yCntr(6,17)
zCO(inactive)=zCntr(6,17)
bCO(freeze)=B(6,17)
xNH(inactive)=XCntr(2,18)
yNH(inactive)=yCntr(2,18)
zNH(inactive)=zCntr(2,18)
bNH(freeze)=B(2,18)
distCOandNH=sqrt[(xCO-xNH)^2+(yCO-yNH)^2+(zCO-zNH)^2]*0.529
distCOandNH(NSteps=30, StepSize=-0.10)

The first line of the GIC input defines the x-coordinate of the centre of the CO bond, and the following two do the same for y and z. They are set to inactive as we only need them as an intermediate anchor for the distance. They will therefore not show up in the output.
The fourth line freezes the CO bond length. I decided to include that because it is part of your initial set-up. After running the calculation I think it should be relaxed though.
The four lines after that define the same as above for the NH bond.
The second to last line defines the distance between the centres. The multiplication with 0.529 converts from bohr to angstrom. That is basically only necessary because it makes the following operation, the actual scan in the last line, easier to specify.

This new interface therefore provides many wonderful opportunities, and possibilities for the most exotic scans. However, it also lets you define certain values to be printed in the output, for much easier post-processing. I can only recommend reading the manual and making a lot of use of this new system.

In order to not leave you hanging dry, here is the resulting scan animated:

animated scan

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you for your answer. What if one wants to perform the equivalent rigid scan ? I could not find a straightforward solution without employing Z-matrix ... $\endgroup$
    – gluuke
    Commented Jul 27, 2020 at 10:45
  • $\begingroup$ @gluuke I don't think there is an easy straight forward way to do a rigid scan without a z-matrix. Of course you can created all resident coordinates yourself and then freeze the majority, but that is probably as extensive a task as creating a sane matrix. $\endgroup$
    – Martin - マーチン ♦
    Commented Jul 27, 2020 at 19:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.